Global genetic diversity, introgression, and evolutionary adaptation of indicine cattle revealed by whole genome sequencing

Domain: earth_science
Published: 2023-11-28
Editorial outcome: accept

Abstract
Indicine cattle, also referred to as zebu (Bos taurus indicus), play a central role in pastoral communities across a wide range of agro-ecosystems, from extremely hot semi arid regions to hot humid tropical regions. However, their adaptive genetic changes following their dispersal into East Asia from the Indian subcontinent have remained poorly documented. Here, we characterize their global genetic diversity using high-quality whole-genome sequencing data from 354 indicine cattle of 57 breeds/populations, including major indicine phylogeographic groups worldwide. We reveal their probable migration into East Asia was along a coastal route rather than inland routes and we detected introgression from other bovine species. Genomic regions carrying morphology-, immune-, and heat-tolerance-related genes underwent divergent selection according to Asian agro-ecologies. We identify distinct sets of loci that contain promising candidate variants for adaptation to hot semi-arid and hot humid tropical ecosystems. Our results indicate that the rapid and successful adaptation of East Asian indicine cattle to hot humid environments was promoted by localized introgression from banteng and/or gaur. Our findings provide insights into the history and environmental adaptation of indicine cattle.

Introduction
The domestication of aurochs (Bos primigenius) gave rise to two distinct but cross-fertile cattle subspecies, humpless taurine cattle (B. taurus taurus) and humped indicine or zebu cattle (B. t. indicus)1. These events were important in human history, with extensive implications for the diet, culture, and socioeconomic structure of the farming populations across the Old World. Taurine cattle were domesticated ~10,000 years before present (YBP), followed by the domestication of indicine cattle 2000 years later in the Indus Valley of modern Pakistan2. Indicine cattle are recognized by their thoracic hump, low metabolic rate, many large sweat glands, large skin surface, and short smooth coat3. They are often resilient to local ticks and capable of tolerating the hot and/or humid climates of the semi arid and tropical regions3. Thus, they can experience a much larger complementary thermal stress spectrum than taurine cattle, which are notably absent in the tropical areas of Asia. The latter are indeed largely confined to temperate to cold environments, with the exception of the West African taurine cattle living in humid and sub humid, tsetse fly-infested, and tropical environments4.

Indicine cattle are the most abundant and important livestock species in South Asia, East Asia, and Africa5, and they represent more than half of all cattle populations worldwide3. The successful global and agro-ecological dispersal of indicine cattle is unique among domestic bovine species. It has been essential for the development of local agricultural lifestyle and economy that have shaped modern societies in subtropical and tropical regions3,4. Their adaptation to the hot climate will be increasingly important in the context of climatic changes, with increasing temperatures affecting livestock production worldwide4.

Archaeological evidence indicate the presence of domesticated indicine cattle earlier in the Indus Valley (~8000 YBP) than in South India (~5000 YBP) and the middle Ganges (~4000 YBP)5,6. The global dispersal of indicine cattle started in the Indus Valley at ~5000 YBP, followed by their spread into Southwest and Central Asia, East Asia, and Africa between 4000 and 1300 YBP5,6. An ancient DNA analysis indicated widespread male-mediated introgression of indicine cattle from the Indus Valley into the Near East from 4200 YBP. Modern DNA analyses have now well documented this male-mediated indicine admixture into African taurine cattle in the eastern, western, and southern areas of the continent5. There is also small but significant indicine introgression into almost all southeastern European cattle breeds7. The expansion of indicine cattle has continued until recent times, and indicine cattle imported in the nineteenth and twentieth centuries into America and Australia have formed large local populations8. Along with their global spread, admixture with local taurine cattle, wild and/or domesticated banteng and gaur, and possibly other unsampled wild bovine species supposedly led to the diversification of indicine cattle populations5,9,10. A common practice is to hybridize other bovine species with cattle to rapidly improve their adaptation to new environments. The establishment of stable hybrid populations is difficult because hybrid males are often sterile, but limited introgression after backcrossing several generations of female hybrids to male cattle is possible11.

Accordingly, three major domestic indicine autosomal lineages are recognized today: (1) the source population in South Asia; (2) African indicine cattle admixed with African taurine diversity5; and (3) East Asian indicine cattle12. Global indicine diversity is further characterized by two Y chromosome haplogroups (Y3A and Y3B)9, two major indicine mtDNA haplogroups (I1 and I2) in Asia6, taurine mtDNA haplogroups in African, American, and Australian indicine cattle populations, and banteng mtDNA in several Indonesian indicine breeds13. Taken together, current autosomal, Y-chromosomal, and mitochondrial ancestries indicate complex domestication and evolutionary processes in the formation of global indicine cattle diversity.

The aim of this study was to explore the unique genomic characteristics and phylogeographic patterns of the diversity of indicine cattle using the largest indicine cattle genome dataset available to date. We present a comprehensive genomic analysis of the variations in the Y chromosome, mitogenomes, and whole nuclear genomes of 354 indigenous indicine cattle sampled from 57 breeds/populations representing the majority of indicine cattle groups. Our findings reveal a discontinuous geographic pattern of genomic diversity and extensive introgression of banteng and gaur, and provide insights into the genomic background of the unusual physiological features that enable indicine cattle to tolerate extreme environments (hot-dry and hot-humid) and a high infectious disease burden.

Results
Genetic diversity and differentiation of indicine cattle
A total of 297 new genomes, including 287 indicine cattle representing 42 breeds/populations and 10 taurine cattle representing three breeds, were sequenced to an average depth of 11.72×. They were combined with 198 (67 indicine and 131 taurine genomes) publicly available genomes (Fig. 1a, Supplementary Note 1, and Supplementary Data 1). Twenty-two whole genomes from other wild and/or domestic bovine species (five gaur, eight banteng, two bison, two wisent, three yak, and two swamp buffaloes) were included for introgression analysis. Sequence reads were aligned to the taurine cattle reference genome (ARS-UCD1.2) and Btau5.0 Y chromosome with an average alignment rate of 99.50% and a coverage of the reference genome of 94.76% (Supplementary Data 1). A total of 354 indicine genomes representing 57 breeds/populations and 141 taurine genomes from 17 breeds/populations were classified as follows: African taurine (AFT, n=19), European taurine (EUT, n=62), Eurasian taurine (EAT, n=28), Tibetan taurine (TBT, n=8), Northeast Asian taurine (NEAT, n=24), African indicine (AFI, n=111), South Asian indicine (SAI, n=118), Southeast Asian indicine (SEAI, n=28), Tibetan indicine (TBI, n=7), Southwest Chinese indicine (SWCI, n=4), East Asian indicine (EAI, n=80), and American indicine (AMI, n=6) cattle (Fig. 1 and Supplementary Data 1). A total of 67,162,108 autosomal SNPs were identified (Supplementary Tables 1 and 2).

The genome-wide nucleotide diversity revealed by autosomal SNPs was generally higher in indicine cattle (0.00261–0.00337) than in taurine cattle (0.00136–0.00164). The highest value (0.00337) was observed within the EAI cattle (Supplementary Fig. 1 and Supplementary Note 2), while the average values for AFI and SAI cattle were 0.00265 and 0.00261, respectively. The level of inbreeding measured by runs of homozygosity (ROH) was lower in indicine cattle than in taurine cattle (Supplementary Fig. 2). Genetic distances estimated via the pairwise fixation index (FST) ranged from 0 to 0.205 between indicine breeds/populations and from 0.201 to 0.550 between indicine and taurine groups (Supplementary Fig. 3).

Principal component analysis (PCA) of the autosomal SNPs revealed clear phylogeographical differentiation, with PC1 corresponding to the contrast between indicine and taurine cattle (Fig. 1b). The PCA and phylogenetic tree almost completely separated the three indicine geographic groups of SAI, AFI, and EAI cattle (Fig. 1c and Supplementary Figs. 4–6). SWCI fell in genetically intermediate positions between SEAI and EAI cattle. The indicine cattle of TBI (Tibet, China) and Nepal were close to SAI cattle (Fig. 1c). The ADMIXTURE analysis recapitulated a similar pattern and identified three indicine and three taurine phylogeographic groups at K=6 (Fig. 1d, Supplementary Fig. 7, and Supplementary Table 3). The same differential topology was observed in the maximum likelihood (ML) tree of these breeds/populations (Supplementary Fig. 8).

Adaptation of indicine cattle
Outside the monsoon season, the Indus Valley has a semiarid climate with a high temperature, high solar radiation, and low rain fall14,15. Because it is the center of origin of indicine cattle, these cattle may be expected to be particularly adapted to such environmental conditions. This has driven the successful spread of indicine cattle into the central and southern regions of the globe. To identify these ancestral adaptations at the genome level, we combined SAI, AFI, and EAI populations for a comparison with taurine cattle by using FST, θπ ratio, and cross-population extended haplotype homozygosity (XP-EHH) approaches (Table 1, Supplementary Note 3, Supplementary Table 4, and Supplementary Figs. 9 and 10). A total of 156 nonoverlapping windows of 50 kilobase (kb) in size were detected using all three approaches. They overlapped with 117 candidate genes (Supplementary Table 4).

The top selection signatures were in two regions on Bos taurus autosome (BTA)7, together spanning 4.46 megabases (Mb) (43.04–44.67 and 50.14–52.97 Mb) (Supplementary Fig. 11). This region was previously identified to be associated with host immunity, environmental thermal stresses, and reproduction in African humped cattle5 (Table 1), supporting its ancestral indicine origin. Another strong selection signature for a gene-rich region was located on BTA19, spanning 1 Mb and covering genes related to antiviral immunity (SPAG716), neurodegenerative disease (KIF1C17), skeletal development (PFN118), cardiac growth (CAMTA219), and muscle development and glycogen storage (ENO320) (Supplementary Fig. 12). We also identified a functional gene, LIPH, on BTA1 (81.58–81.69 Mb) (Supplementary Fig. 13), which was linked to hair growth deficiency in humans21, implying its potential contribution to the heat tolerance of indicine cattle via the control of coat hair length and/or thickness. Seven of the remaining 75 genes in the topmost significant sweep regions are functionally associated with heart development, blood circulation, DNA damage, and light response (Table 1 and Supplementary Table 4). However, further research is warranted to test their roles in heat adaptation or other differences between indicine and taurine cattle.

Adaptation of South Asian indicine cattle
Throughout the history of migration and admixture of indicine cattle, genomic regions under selection might have been lost in specific indicine groups. We therefore also performed a test for positive selection signatures in SAI cattle using θπ, iHS, CLR, and FST estimates based on the comparison between SAI and non-SAI groups (Supplementary Data 2). SAI comprised ancestral indicine cattle that have adapted to the harsh conditions with high temperatures and solar radiation but low rain fall. Gene Ontology (GO) category and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses revealed major enrichments of positively selected genes (PSGs) involved in environmental adaptation, including arachidonic acid metabolic process (GO:0019369) and infectious disease (human papillomavirus infection, bta05165) (Supplementary Table 5). Notably, the arachidonic acid metabolic pathway has been considered as an important molecular mechanism for desert adaptation22,23 via efficient water reabsorption in sheep24, chickens23, and Bactrian camels22, implying a deep convergent evolution among indicine cattle and other species survival over high salt diets in similarly arid and hot environments. We also detected several PSGs in SAI cattle, to be associated with fluid homeostasis (APELA25) and sensory ability (CALB226) (Supplementary Data 2). Therefore, this comparison revealed a variety of important genes, pathways, and GO categories associated with the genetic adaptation of SAI cattle to semiarid environments (Supplementary Table 5).

Additionally, we compared the light- and dark-coated SAI breeds, such as the white-coated Bhagnari and Dajal cattle from Pakistan, by using FST and θπ ratio estimates. We identified shared selective sweeps around pigmentation loci, e.g., LEF1 and ASIP, in the light-coated indicine breeds (Fig. 2). This selection pressure may have been favored or driven by high temperatures and intense solar radiation and/or human preferences. Across the whole genomes, the CLR and iHS analyses revealed 368 regions overlapping with 477 genes present in AFI and SAI cattle (Supplementary Data 3), supporting that the ancestral adaptations of SAI cattle were equally important for AFI cattle.

The indicine adaptation to the tropical, humid Asian environment
Although domestic indicine cattle emerged in South Asia6, different phylogeographic groups have evolved to adapt to local environments during their global dispersal, including nonnative humid, tropical regions27. The agroecologies of southern East Asia are referred to as mixed subhumid or humid systems in contrast to the mixed arid or semi arid system of South Asia28. Southern East Asian agroecologies are characterized by high humidity and rainfall, as well as a relatively high incidence of tropical diseases10,28.

We first assessed whether any ancestral indicine genomic regions under selection were present in EAI cattle. We used the population branch statistic (PBS) with banteng as an outgroup to detect recent selection signatures in EAI cattle, while avoiding the effect of introgression from other bovine species on the selection of EAI cattle (Fig. 3 and Supplementary Table 6).

The longest gene-rich region of 410 kb under selection was observed on BTA22 (Fig. 3), which contained 14 protein-coding genes. Two of these genes are related to the host immune system (MST1R29,30 and MON1A31) and four tumor suppressor genes are related to lung cancer (SEMA3B32, GNAI233, SEMA3F34, and RBM535) (Fig. 3). Moreover, TRPA1 on BTA14 was highly differentiated in EAI cattle and may be associated with thermal avoidance behavior (Fig. 3b). Other PSG associated with circannual clock (FAM204A36) may help EAI cattle adapt to seasonal changes in day length by modifying their behavior such as grazing at night. In addition, PSGs involved in the immune system (Supplementary Table 6) could confer tolerance to parasitic pathogens in EAI cattle.

We then assessed whether the local adaptation of EAI cattle may have involved introgression from other Asian bovine species. Increasing evidence for past introgressions between bovine species distributed in East Asia and Southeast Asia is available9,10,12,37. Such introgression may have facilitated the rapid adaptation of indicine cattle to humid regions9. Previously, we reported the introgression from banteng into southern Chinese cattle9, while another study suggested admixture from an unknown source into EAI cattle12. In the present study, we used TreeMix38, RFmix39, the D statistic40, and phylogenetic analyses of candidate introgressed fragments to investigate the gene flow from banteng and gaur into 97 EAI and SWCI genomes (Supplementary Note 4, Supplementary Table 7, Supplementary Figs. 14–16, and Supplementary Data 4 and 5). The proportions of banteng and gaur ancestries ranged from 1.13% to 10.21% and from 2.06% to 9.98% in the EAI genomes (Fig. 4a and Supplementary Data 6 and 7). EAI cattle in the southeastern coastal region of China showed the highest level of banteng and gaur ancestries (Supplementary Fig. 17). We used the U20 statistic to identify frequently introgressed genes in the EAI genomes41. We calculated U20 (1%, 20%, and 100%) to be equal to the number of alleles where banteng or gaur had a particular allele fixed, while its frequency was less than 1% in SAI but greater than 20% in EAI genomes (Supplementary Fig. 18)41. We found 1267 genes in the EAI genomes to be of banteng origin and 1488 genes to be of gaur origin, with 921 genes shared by both banteng and gaur (Supplementary Fig. 19). GO analysis revealed significant overrepresentations of introgressed genes involved in biological processes contributing to environmental adaptation, the nervous system, and the endocrine system (Fig. 4b and Supplementary Table 8).

Using a higher cutoff for the frequencies of banteng- or gaur-derived alleles in the EAI genomes (U50 (1%, 50%, 100%)), 70 introgressed genes in 32 candidate regions were shortlisted (Fig. 4c and Supplementary Tables 9 and 10) and then validated by phylogenetic analyses (Supplementary Figs. 20–25). Notably, one region on BTA1 (66.70–66.80 Mb) harbored genes relevant to water transport (ILDR1)42, blood homeostasis (CASR)43, and skin disease (CSTA)44 (Fig. 4c and Supplementary Fig. 26). Another region on BTA25 (0.21–0.26 Mb) also demonstrated a clear pattern of introgression in EAI cattle (Fig. 4d–g). The geographic distribution of haplotypes in global cattle populations showed that the introgressed haplotype had the highest frequency in EAI cattle, which was supported by phylogenetic analysis (Fig. 4d), FST values (Fig. 4e), and the degree to which the haplotypes were shared (Fig. 4f). This region contained a cluster of genes (HBM, HBA, HBA1, and HBQ1) involved in biologically relevant oxygen transport (GO:0015671) (Fig. 4 and Supplementary Table 8), which were also associated with resistance to severe malaria in humans45. Within the genes in the hemoglobin family cluster, 11 missense mutations showed significantly altered frequencies of specific alleles between the EAI and other indicine groups (Fig. 4g, h).

Uniparental dispersal of indicine cattle
Specific uniparental lineages may be informative for the reconstruction of historical migration patterns. Here, we identified 1389 SNPs in the male-specific region of the bovine Y chromosome (MSY) in 309 males (Supplementary Data 1), which were defined as the taurine haplogroups of Y1, Y2A, and Y2B and indicine haplogroups of Y3A and Y3B reported previously9 (Supplementary Note 5 and Supplementary Figs. 27–30). Within Y3A, we resolved two minor sub-haplogroups of Y3A1 and Y3A2 and a major sub-haplogroup of Y3A3, whereas Y3B was divided into the sub-haplogroups of Y3B1, Y3B2, Y3B3, and Y3B4 (Fig. 5a). Most of these sub-haplogroups were present in SAI cattle, supporting that South Asia was a primary center of paternal genetic diversity of indicine cattle. Following a west–to–east genetic cline, the haplogroup Y3A was predominant in EAI (88.13%) and North-Central Chinese cattle (Supplementary Fig. 29)9. In contrast, Y3B was dominant in SAI and AFI cattle (89.76%) (Fig. 5a). A phylogenetic tree showed the divergence of Y3A1, followed by Y3A2 and then Y3A3, which correlated with their geographic ranges: Y3A1 in SAI and SEAI; Y3A2 in TBI and SWCI; and Y3A3 as the predominant sub-haplogroup in EAI and North-Central Chinese cattle (Supplementary Data 1 and 8 and Supplementary Fig. 29). Y3A3, which occurred only in Indochina and China, may have emerged as a new sub-haplogroup during the indicine eastward migration. Among the Y3B haplogroups, Y3B2 migrated to the east, while Y3B4 had a larger range and was almost exclusively present in AFI cattle. An ML tree of 347 assembled mitogenomes separated taurine and indicine cattle first and then added indicine I3 as a new haplogroup and I1a as a recent split-off of I1 (Fig. 5b and Supplementary Figs. 31–35). Similar to the Y3A3 sub-haplogroup, I1a may have emerged as a new sub-haplogroup during the indicine eastward migration and got established in EAI and some North-Central Chinese cattle (Supplementary Data 8 and Supplementary Fig. 34)9. Notably, both the westward and eastward migrations led to the fixations or near fixations of unique indicine Y-chromosomal and mitogenome haplogroups (Fig. 5a, b).

Demographic history of indicine cattle
We used the multiple sequential coalescent Markovian model (MSMC) to infer historical changes in the effective population size and population separation of the three core indicine groups (EAI, SAI, and AFI). All indicine groups diverged from taurine groups between 251.5 and 301.2 thousand years ago (kya) and experienced common and substantial declines in Ne at 20–30 kya and 7–9 kya. We observed earlier and clear splits at ~10.3 kya between EAI and SAI cattle and at ~11.8 kya between AFI and SAI cattle (Supplementary Fig. 36).

Using an empirical Bayesian approach, we detected more recent estimates of divergence times among indicine Y haplogroups than for mitogenome lineages. The divergence of Y3A from Y3B occurred at ~23.1 kya (Supplementary Fig. 30), while the newly emerged Y3A3 diverged from both Y3A1 and Y3A2 at ~20 kya (Supplementary Fig. 30). The Bayesian skyline plot (BSP) of indicine Y chromosome haplogroups displayed a slow population expansion from 10 kya, followed by an evident increase from 5 kya, which overlapped with the expansion of the post-domestication indicine populations (Supplementary Fig. 30). The divergence of sub-haplogroups within indicine mitogenomes occurred from 9.6 to 24.8 kya (Supplementary Fig. 32). The BSP of indicine mitogenomes showed a population decrease after domestication (10–8 kya) but a rapid increase from 7 kya, which overlapped with the expansion of the indicine populations after domestication (Supplementary Fig. 35). We also observed a population increase for I1a at ~3.73 kya, due probably to the expansion of I1a into the current distribution in East Asia (Supplementary Fig. 35).

The Southeast Asian coast was the main entry point of indicine cattle into East Asia
Indicine cattle were believed to have entered East Asia from South Asia through the inland routes8,9. However, we did not observe any west-to-east genetic clines, but a rather abrupt transition, as evidenced by different genetic features: (1) the drastic shift of haplogroup frequencies of both uniparental markers, with the Y-chromosomal sub-haplogroup Y3A3 and the maternal sub-haplogroup I1a exclusively found in the Indochinese and Chinese cattle populations (Fig. 5); (2) the autosomal variation indicated relatively long-distance dispersals between SEAI and SWCI and between SWCI and EAI cattle (Fig. 1c); and (3) the extent of wild and/or domestic bovine introgression into EAI genomes did not follow a gradual west-to-east cline (Supplementary Fig. 17). Remarkably, all these transitions were geographically close to each other across Southwest China, where the mountainous landscape is traversed by three major rivers flowing from north to south (Nujiang River, Honghe River, and Lancang River; Supplementary Fig. 37). We propose that during the indicine eastward migrations, these geographic barriers were circumvented by maritime migrations along the coast (Fig. 5c), as maritime migrations also played an important role in the arrival of indicine cattle in Africa5.

Methods
Ethics statement
Blood samples and ear tissues were collected during routine veterinary treatments with the logistical support and agreement of relevant agricultural institutions in each country. All procedures involving sample collection and experiments were approved by the Animal Ethical and Welfare Committee, Northwest A&F University (Approval No. DK2022065).

Read mapping and SNP calling
We generated genotype data following the 1000 Bull Genomes Project Run8 guideline (http://www.1000bullgenomes.com/) (Supplementary Note 1). We removed low-quality bases and artifact sequences using Trimmomatic v0.3952, and all clean reads were mapped to the taurine cattle reference assembly (ARS-UCD1.2) and Btau_5.0.1 Y using BWA-MEM v0.7.13-r1126 with default parameters53. We then used SAMtools v1.954 to sort bam files. For the mapped reads, potential PCR duplicates were identified using “MarkDuplicates” of Picard v2.20.2 (http://broadinstitute.github.io/picard). “BaseRecalibrator” and “PrintReads” of the Genome Analysis Toolkit (GATK, v.3.8-1-0-gf15c1c3ef)55 were used to perform base quality score recalibration (BQSR) with the known variant file (ARS1.2PlusY_BQSR_v3.vcf.gz) provided by the 1000 Bull Genomes Project.

For SNP calling, we created GVCF files using “HaplotypeCaller” in GATK with the “-ERC GVCF” option. We called and selected candidate SNPs from these combined GVCF files using “GenotypeGVCFs” and “SelectVariants”, respectively. To avoid possible false-positive calls, we used VariantFiltration of GATK as recommended by GATK best practices: (1) SNP clusters with “-clusterSize 3” and “-clusterWindowSize 10” options; (2) SNPs with mean depth (for all samples) < 1/3× and > 3× (×, overall mean sequencing depth across all SNPs); (3) quality by depth, QD < 2; (4) phred-scaled variant quality score, QUAL < 30; (5) strand odds ratio, SOR > 3; (6) Fisher strand, FS > 60; (7) mapping quality, MQ < 40; (8) mapping quality rank sum test, MQRankSum < −12.5; and (9) read position rank sum test, ReadPosRankSum < −8 were filtered. We then filtered out non biallelic SNPs and SNPs with missing genotype rates > 0.1. Imputation and phasing of SNPs were simultaneously performed using BEAGLE v4.0 with default parameters, and SNPs were filtered with DR2 < 0.956. The remaining SNPs were annotated according to their positions using SnpEff v4.357.

Genetic diversity and population genetic structure
The genome-wide nucleotide diversity of different cattle geographic groups was estimated with VCFtools v0.1.1658. Genetic distances between breeds/populations were calculated with the FST estimates and ROH were analyzed using PLINK v1.959,60 (Supplementary Note 2). For PCA and admixture analysis, we first filtered out SNPs with a minor allele frequency (MAF) < 0.01 and performed LD-based pruning for the genotype data using the --indep-pairwise 50 10 0.1 option of PLINK v1.957. For PCA, we used the Smartpca program in EIGENSOFT v4.261. The Tracy-Widom test was used to determine the significance level of the eigenvectors. ADMIXTURE v1.3.0 was used to quantify genome-wide admixture among cattle breeds/populations62 and run for each possible group number (K=2 to 8), where K was the assumed number of ancestries. The delta K method was used to choose the optimal K62. A neighbor-joining (NJ) tree was constructed using the matrix of pairwise genetic distances calculated by PLINK v1.957. A population-level phylogeny was reconstructed using the maximum likelihood (ML) in TreeMix38.

Detection of selection signatures shared by all indicine cattle groups
We screened genomic regions under selection with the largest differences in genetic diversity (θπ (indicine/taurine) ratio) and FST outliers between taurine (EUT, EUAT, NEAT, TBT, and AFT, n=141) and indicine (SAI, EAI, and AFI, n=309) cattle using VCFtools v0.1.1658 (Supplementary Note 3). We also performed XP-EHH analysis using the default settings of selscan v1.163. The π ratio, FST, and average normalized XP-EHH score were calculated for 50 kb windows with 20 kb steps. Top 1% windows were identified as significant genomic regions.

Detection of selection signatures in South Asian, East Asian, and African indicine cattle
The CLR and iHS were employed to detect these selection signatures in the SAI (n=118), EAI (n=80), and AFI (n=111) genomes. The CLR was calculated for sites based on 50 kb windows with 20 kb steps using the SweepFinder264 command “SweepFinder2 -lu GridFile FreqFile SpectFile OutFile”. The iHS was calculated in selscan v1.163, and the proportion of SNPs with |iHS| ≥ 2 was calculated in windows of 50 kb and steps of 20 kb. To perform iHS and CLR computation, information on the ancestral and derived allele states is needed for each SNP. In our analysis, the ancestral allele was defined as the allele fixed in the swamp buffalo that was included in the genotype call set, while the SNPs failed in genotyping call for their ancestral state were discarded. To capture potential genes that were specifically selected for each indicine group, we also calculated the FST between the target group and two other indicine groups. p values were calculated for the CLR, |iHS|, and FST windows, and the overlap windows of p < 0.005 (Z test) of each method were considered candidate signatures of selection.

Considering that the EAI genomes were affected by banteng/gaur introgression, we used the PBS65 in 50 kb windows with 20 kb steps to scan for genomic regions highly differentiated in EAI relative to SAI, AFI, and banteng (n=4) genomes. Significant genomic regions were identified by a p < 0.005. In addition, FST and θπ methods were used to visualize the line chart of the top signals.

Introgression analysis
TreeMix38, the D statistic40, and RFMix v2.0239 were used to determine the gene flow between EAI and other bovine species. OptM was used to determine the optimal number of migration edges in the TreeMix. RFMix was used to identify regions introgressed from banteng or gaur into EAI cattle39 (Supplementary Note 4). Pure taurine cattle, SAI cattle, banteng or gaur were selected as the reference panel according to D and f3 statistics (Supplementary Note 4). We calculated the probability of banteng/gaur introgressed tracts in EAI cattle due to incomplete lineage sorting (ILS)66. We let r be the recombination rate per generation per base pair (bp) in indicine cattle, m be the length of the introgressed tracts, and t be the length of other bovine species (banteng and gaur) and cattle branches since divergence10. The expected length of a shared ancestral sequence was L = 1/(r × t) = 206.52 bp. The probability of a length of at least m was 1-GammaCDF(m, shape=2, r=1/L), in which GammaCDF is the gamma distribution function. We applied the probability of ILS < 0.05 to filter short introgressed segments in the RFMix results (Supplementary Data 6 and 7). A total of 79 topological trees were used to confirm banteng or gaur introgression and visualized by DensiTree67 (Supplementary Fig. 16). Second, we used the statistics U20 (1%, 20%, and 100%) and U50 (1%; 50%; 100%) to detect sites based on 50 kb windows with 20 kb steps where banteng or gaur had a particular allele at a frequency of 100%, while the frequency was less than 1% in SAI but greater than 20% or 50% in EAI cattle (Supplementary Tables 9 and 10).

Functional enrichment analyses
As source of annotation, we used the source Bos taurus Annotation Release 106 (GCF_002263795.1_ARS-UCD1.2_genomic.gtf) based on the NCBI assembly of GCF_002263795.1. Gene set enrichment analyses were carried out with GO categories and KEGG pathways for KOBAS v3.068. The value was calculated using a hypergeometric distribution and corrected for the FDR. To adjust for multiple testing, pathways with p < 0.05 were considered significantly enriched.

Paternal analysis
The X-degenerate region that consists of single-copy genes within the male-specific part of the Btau_5.0.1 Y chromosome reference sequence (GCF_000003205.7) was selected (Supplementary Note 5). After removing sites with missing genotypes in 10% of the samples, 1389 SNPs were extracted. Fasta sequence files were used to generate haplogroup trees. Phylogenetic construction was performed using BEAST v2.6.069. To further explore the migration of Y3A haplotypes in China, we extracted 26 indicine Y haplotypes representing 11 hybrid breeds from North-Central China in previous studies9. We genotyped 26 individuals according to these 1389 SNPs (Supplementary Data 8).

Bayesian age estimates of haplogroups and Bayesian skyline plots (BSPs) were generated using BEAST v2.6.070. A maximum clade credibility tree was generated using a 10% burn-in with TreeAnnotator and drawn with FigTree v1.4.371. The BSPs of the indicine sub-haplogroup Y3 and its sub-haplogroup Y3A3 were generated. The following parameters were applied in both runs: HKY substitution model with gamma-distributed rates, a log-normal relaxed clock, coalescent Bayesian skyline analysis, a mutation rate per generation of 1.26 × 10−8, and a generation time of 6 years72. The node age of sub-haplogroup Y3A3 (5.57 kya) was used as the only a priori parameter. We ran 100,000,000 iterations for Y3 and 50,000,000 iterations for Y3A3, with samples collected every 5000 steps, and visualized the BSPs obtained with Tracer v1.773. LogCombiner was used to perform 10% burn-in. The results were calibrated with a generation time of 6 years, and BSP plots were plotted using the ggplot2 in R v4.1.074.

Mitogenome phylogeny
We assembled and selected 344 mitogenomes and aligned them to 18 bovine mitogenomes (Supplementary Note 5). Phylogenetic relationships were inferred from the mt genomes using RAxML v8.2.975 with the following parameters: -f a -x 123 -p 23 -# 100 -k -m GTRGAMMA. The final tree topology was visualized using FigTree v1.4.371. The median-joining network was constructed using NETWORK v5.0.1.171. We extracted mitogenomes representing 13 hybrid breeds from North-Central China to further explore the migration of I1a sub-haplogroup in East Asia9 (Supplementary Data 8).

Bayesian age estimates of haplogroups and BSPs were generated using BEAST v2.6.070. BEAST runs were performed on three datasets with mitogenome coding regions (all 362 mitogenomes, 243 indicine mitogenomes, and 119 taurine mitogenomes). We used the HKY substitution model (with gamma-distributed rates) with a log-normal relaxed clock. We applied an evolutionary rate of 2.043 ± 0.099 × 10−8 base substitutions per nucleotide per year76. For each dataset, we performed ten independent BEAST runs with the chain length established at 20,000,000 iterations, samples collected at every 5000 MCMC steps and applying a 10% burn-in. The runs were then combined using the LogCombiner utility within the BEAST package by applying another 10% burn-in. A maximum clade credibility tree was drawn with FigTree v1.4.371. BSP data were obtained with Tracer v1.7.173 using default parameters and calibrated using a generation time of 6 years77. The BSPs were plotted using the ggplot2 in R v4.1.074.

Estimation of effective population size and divergence time using autosomal SNPs
The multiple sequential coalescent Markovian model 2 (MSMC2)78 method was used to model the population history of the three core indicine groups (EAI, SAI, and AFI) and to infer historical changes in their effective population size and population separation (Supplementary Note 5).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "General comments\nChen et al. present the results of whole-genome resequencing of 355 indicine cattle\ngenomes and 141 taurine cattle genomes from 57 and 17 populations, respectively. The\nquantity and quality of the material studied here, the methods used and the results\nobtained are of a high standard and deserve publication in a high-ranking journal. However,\nthe main problem concerns the core zebu breed group and the outgroups.\nRegarding the core Zebu group: The Zebu cattle breeds from India were missing in this\nstudy. The surrogates used for this purpose are Gir, Nellore and Brahman, which were bred\noutside India for centuries in new environments and with sporadic interchange with Bos\ntaurus breeds. Introgression of Bos taurus into the WGS of Gir, Nellore and Brahman used\nhere is neither ruled out nor investigated by the authors. Brahman in particular is known to\nbe a synthetic breed with Bos taurus content. Many countries such as Sri Lanka, Pakisitan\nand Nepal (SAI group) use Bos taurus breeds to improve local cattle breeds. Therefore, the\npurity of any reference used for introgression studies should be tested first.\nWith regard to outgroups: Besides bison and yak, two Bantengs and two Baurs are used for\nimportant introgression analyses in this study. Again, there is a possibility of introgression of\ncattle into Banteng and/or Baurs. This possibility has not been ruled out and has not been\ninvestigated. There are some unpublished results confirming cattle DNA segments in\nBantengs. I do not have comparable information for Gaur, but this is possible. As shown in\nthe study analysing cattle introgression in Yak, RFMix can produce false positive signals if\nthe reference population is itself introgressed.\nThe next problem is the sample size of the outgroups (Banteng and Guar) in the\nintrogression analysis. I expect that the relatively small sample size of the outgroups (2\nsamples each) will affect the performance of RFMix and the determination of ancestral\nalleles and their frequencies in the U20 and U50 analyses. In fact, the original paper of\nRFMix already describes how the sample size of the reference panel affects its performance\n(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3738819/figure/fig6/). How did the\nauthors determine the robustness of their results?\nThe results of the D-statistics are confounded by the structure of the ancestral population as\nwell as incomplete lineage sorting. I recommend that the authors perform other\nindependent approaches, e.g. those based solely on branch length and tree discordance, to\nvalidate these results.\nHow did the authors determine the optimal number of migration edges in the TreeMix?\nIdeally, they should also include the matrix of residual values to show how integrating\nmigration edges improved the % of variance explained.\nHow did the authors determine the expected length of introgressed haplotypes? How did\nthey make sure that it is significantly different from the expected length because of its\nshared ancestral origin? Does the length of introgressed haplotypes tell us anything about\nthe age of introgression?\nEach figure and table, including supplementary figures and tables, should be considered as a\nstand-alone object, i.e. the readers should be able to understand it without reading other\nobjects. For this reason, the legend of the figure or table must be informative enough to\nexplain the main message of this stand-alone object. This is not the case. I make a specific\ncomment for one figure only (mentioned in the comments below), but it applies to all and is\ntherefore a general comment.\nSpecific comments\nLine 87-88, “..resilient to local pathogens”: consider replacing with tick resistance as also\nmentioned in the reference paper cited here.\nLine 106-108: Please keep in mind that there is a small but significant Bos indicus\nintrogression also into almost all south-eastern European cattle breeds up to the southern\nAlps (https://doi.org/10.1186/s12711-020-00560-8).\nLine 108-109, “….Gujarat in Western India”: While the reference paper do discuss the\nhypothesis about possible wild introgression in southern Indian zebu, the same is not\nmentioned for Gujarat, consider rewriting this sentence.\nLine 144-145, “thirteen whole genome sequencing data…”: The authors should mention\nhere the number of WGS used for each species separately.\nLine 166-167, ”…the PCA and a phylogenetic tree..”: The authors should describe here how\nSWCI fits into the overall phylogenetic reconstruction.\nLine 176-178, “We identified candidate selected genomic regions…”: Why author only\ncombined the three indicine groups (SAI, AFI, EAI) and left the other two (SEAI and SWCI) for\nthis analysis?\nLine 175ff: The authors use reverse genomics to detect candidate genes. They constructed\ntwo groups of animals/populations: 1) taurine and 2) SAI, AFI and EAI indicine populations.\nBoth groups lived in environments that differed by a variety of factors, including\ntemperature. The authors found numerous significantly differentiated regions and attribute\nthem all to adaptation to high temperatures. The idea of reverse genomics is to find the\noutliers and then look for possible causes. Only after some causality analysis can we\nconclude that some of the outliers are an adaptation to high temperatures or something\nelse. The authors should consider significant selection signatures as anonymous candidates.\nThese anonymous candidates are subject to positive selection into a group of animals that\nshare a specific environment. These anonymous candidates could, but need not, be caused\nby one specific environmental factor.\nLine 182-183, “In addition, the expression level…”: As the tissue and animals (with respect\nto gender and age) are quite heterogenous, I will take this result with pinch of salt,\nmoreover, how many other genes were differentially expressed and not under positive\nselection?\nLine 192-195: Again, the authors link LIPH to a possible contribution to heat tolerance in\nindicine cattle by controlling hair length and/or thickness. However, we (readers) are not\ninformed whether indicine and taurine cattle differ in hair length and thickness? I cannot\nfind any such information in this manuscript.\nLine 202ff: The SAI cattle group consists of breeds from Sri Lanka, Pakistan and Nepal (a very\nbroad area), but all outliers are defined as candidates for adaptation to semi-arid\nenvironmental stress.\nLine 223, “The genome-wide analysis shows…SAI origin”: This is well-established fact;\ntherefore, I would remove it or declare as confirmation.\nLine 319-321: With regard to mtDNA distribution, a comparable situation can be observed in\nAFI and EAI. However, the presence of haplogroup T in Africa was not interpreted as T\nhaving emerged as a new haplogroup during the indicine westward migration, but as a\nhaplogroup that was already present in the area before the indicine introgression. The\nauthors should consider both the time of coalescence and the time of migration when\ndrawing their conclusions.\nLine 346-361: Samples originating from India and sampled there are missing in this study.\nThis could be a reason why the west-east cline is missing and an abrupt transition is\nobserved. I did not understand this section of the text. There are several problems in the\ntext itself and in the corresponding figures. The legends of the figures are generally not\ninformative enough. Let's just take Supplementary Figure 17 as an example. The reader\ncannot understand this figure from its legend. Even more, there are many trees labelled\nwith the same A, B, C,.... and we do not know which tree represents which haplotype, etc. I\nalso do not know what to conclude from figure 50, i.e. the legend does not help me to\nunderstand the message from line 355-357.\nLine 414: (in brackets)???\nLine 447-458: Two recent and comprehensive studies on paternal and maternal haplogroups\nin cattle could be helpful here (DOI: 10.1111/eva.13315 and doi: 10.1111/age.13104).\nLine 459-470: This can be understood to mean that East Asia was a land without cattle, so\nthat Bos indicus spread with the single matriline I1a in this empty area. This was also the\ncase with the spread of domestic sheep in Europe. However, for sheep, there is\narchaeological evidence for a sheep-less Europe before the Neolithic. Do you have any\nevidence that East Asia was a cattle-less region before the spread of Bos indicus? Africa was\nnot cattle-less either, and as a result of paternal introgression we only observe haplogroup\nT. Please clarify.\nLine 470-475: The reason behind this could be lack of sampled population along the inland\ntrading routes from South to East Asia. Therefore, the integrated genetic analyses as\nmentioned in the paper may have difficulty supporting inland trading even if this is correct.\nLine 558-559, “Treemix and Dsuite…”: Dsuite is a tool and not any method, please clarify\nthis.\nSupplementary table 1:\nIf possible, the country of origin should also be mentioned, also verify that breeds with\nName “SriLanka” exists.\nSupplementary table 2:\nHariana is declared as an Indian breed and the sample used here was sequenced in Chen et\nal. 2018 and should be from India. We know that India is very restrictive when it comes to\nforeign use of biological resources. There is no Indian collaborator in Chen et al. 2018, so I\ndo not see any legal way to analyze Indian biological material outside India or without an\nIndian collaborator. Was this Hariana cattle sampled in India or elsewhere? Please clarify\nand correct if necessary in Supplementary Table 2.", "reviewer_2": "This manuscript is interesting and provides data from a unique sample of individuals that\nhad not been previously characterized. Overall, I find the manuscript well written, although\nit appears to have been done in parts by different groups of individuals with different\nwriting styles and slightly different nomenclature. This is inevitable with a project of this\nsize. I would encourage the corresponding authors to go through the entire manuscript to\nmake it more uniform. I do have a number of concerns related to the details that were\npresented, or not presented, that make a final determination not possible at this time.\nI always start my reviews with the materials and methods but since there was such a\nsubstantial amount of information in the supplemental information, I will start there.\nUnfortunately, the supplemental document does not have line nor page numbers, so I’ll\nquote specific lines.\nWhole-genome resequencing\n“Additional detailed information on the mapping rate and sequencing depth are provided in\nSupplementary Table 1” In Sup. Table 1 you use the SRR* ID as the sample. This is incorrect\nbecause these are run IDs. You should use the BioSample ID. I downloaded this table and\ncompared these SRR to the NCBI database and you have two duplicate samples.\nSRR2016752 and SRR2016754 are the same SAMN3387026 individual. With only a single\nduplicate, it is unlikely that any of the results would have been significantly affected.\nHowever, this table needs to be corrected to present the BioSample and if data reanalysis is\nperformed, one of these two should be removed.\nPopulation genomic analyses Genetic diversity\n“Runs of homozygosity (ROH) were identified using the --homozyg option implemented in\nPLINK v.1.9, which slides a window of 50 SNPs (--homozyg-window-snp 50) across the\ngenome.” This is a very small window considering the density of markers. It is well\nrecognized that results are sensitive to settings for this analysis (see\nhttps://doi.org/10.1186/s12864-020-6463-x). I question whether or not the settings used\nwere appropriate. Specifically, by using too small of a window size and other inappropriate\nsettings, far too many ROH are detected. Evidence of this can be seen in Supplementary Fig.\n2 where there are a significant number of samples that have estimated ROH totaling greater\nthan 1/3 of the genome and a large number of samples with ROH greater than half the\ngenome. This is unrealistic. This analysis needs to be redone with appropriate settings,\nwhich may need to be determined empirically given the data. It would also be useful at this\ntime to examine any correlation between average genome coverage and number/length of\nROH as one would predict that lower coverage samples will have heterozygous sites\nundercalled which may manifest as ROH when in fact there are not ROH if heterozygous\npositions were called accurately.\nPrincipal component analysis (PCA)\n“We removed \"all LD\" using…” First, please rephrase this. You did not remove “all LD”. I\nhave the same objection here to the parameters used. First, this is a very small window and\nprobably left a significant number of loci with r^2 > 0.20 which are more than 50kb apart.\nGiven everything that we know regarding LD in bovine genomes, you should empirically find\nthe appropriate setting for this analysis to achieve the desired thinning. This does not have\nto be done genome-wide and one could simply use 10Mb from any autosome to estimate\nwhat the appropriate parameters should be. My recommendation would be to use chr25\nsince we know there are no assembly issues. For a 10Mb region, calculate *all* the pairwise\nLD values using ‘-r2 yes really’ and examine the distribution of these values. This figure can\nbe included in the supplement to provide justification for the settings chosen. As the\nmanuscript is currently written, the settings are arbitrary and likely impacted the results and\ninterpretation.\n“PCA was performed on the genome-wide unlinked SNP dataset…” Please change the word\n“unlinked” throughout the manuscript. You do not have any unlinked variants. Perhaps use\n“LD pruned” as a replacement.\nStructure analysis\n“…200 bootstrap replicates were performed…” Why 200? Why not 10, or 10,000. Is the\ndefault of 200 appropriate for your data?\n“…used to determine the optimal ancestry number…” Exactly how did you determine the\noptimal number? There is quite a bit of literature on this, but nothing was stated or cited.\nNeighbor-joining (NJ) and maximum likelihood (ML) phylogenetic trees\n“…A window size of 1,000 SNPs was used to account for linkage disequilibrium…” This\nshould be unnecessary if you properly LD pruned in the first place. This brings up another\nissue I have with the manuscript. You have many different analyses and it is difficult to tell\nexactly which set of loci you are using for which analysis. In this section, it appears that you\nare using the “raw” 67M variants. However, it seems to me that if the objective is to\n“identify closely related individuals” then the LD pruned dataset from the PCA section would\nbe appropriate. However, it appears that these two analyses used different filtering.\nMy recommendation would be to provide a supplemental table that lists the major analyses\nperformed and exactly which data manipulations were done for each analysis. Since data\nmanipulations may impact one analysis but not another, it would be helpful to the reader to\nunderstand exactly what was done and which analyses used the same loci. Furthermore,\nyou can specify the total number of loci in this table for each dataset.\nDetection of selection signatures shared by all indicine populations\n“…were calculated for 50-kb windows with a 20 kb step across the autosomes…” Line 907 of\nthe main manuscript says that a 100kb window was used, which is it? I have issues with\nwhat is specified here. Typically, these types of analyses are performed where the window\nsize is an integer multiple of the step size. For instance 25 kb step and 50 kb window or 20kb\nstep and 100kb window. This needs to be addressed throughout the manuscript and I would\nrecommend using the same parameters for all analyses that use this type of setting (where\nappropriate). This will make comparisons between analyses much easier.\n“…these windows harbored 117 candidate genes in indicine cattle…” I don’t think these 117\ngenes are specific to indicine cattle and you can probably just say cattle. By saying indicine\ncattle here one might think you are implying that these are not in taurine cattle, which isn’t\nthe case.\n“…the expression levels of these 117 candidate genes were validated in nine different\ntissues of taurine and indicine cattle.” I’m not sure exactly what you are validating? The fact\nthat there are RNA-seq data in those tissues for these annotated genes simply validates that\nthe annotation is correct.\n“Candidate genes under selection were defined as those overlapped by sweep regions or\nwithin 20 kb of the signals.” 20kb seems arbitrary. What justification do you have for\nchoosing 20kb? Why not 50kb, or 10 kb or simply require the gene to fall within the sweep\nregion and not allow any overlap. It seems to me that allowing for a gene to be outside your\nsweep region, by some arbitrary distance, speaks to your confidence in defining the\nboundaries of sweep regions. Some justification should be given for whatever value you\nchoose.\n“…were plotted using a 10-kb sliding window…” Again, arbitrary. Why not use the same\nwindow size for plotting that you used for detection?\nDetection of selection signatures in SAI, EAI, and AFI lineages\n“…windows of 50 kb and a step of 20 kb…” Same comment on windows.\n“…the overlap windows of P values less than 0.005 (Z test) of each method were considered\ncandidate signatures of selection…” Why 0.005? Again, some justification is warranted. It\nseems to me that you have an extraordinarily large number of windows that were tested\nwhich implies that some adjustment needs to be considered for multiple testing.\nAppropriately accounting for multiple testing here could significantly change the detected\nregions which in turn could significantly impact all downstream inference based on these\nresults. As written, it is impossible to determine exactly what was done and what impact it\nmay have on the results or interpretation.\n“Only pathways or annotations with a Bonferroni-corrected P < 10−2 were retained\n(Supplementary Tables 7-10).” This should probably be written as P<0.01 to be consistent\nthroughout the manuscript.\nIntrogression analysis\n“Phylogenetic analyses of these segments confirmed the banteng or gaur introgression into\nspecific EAI individuals.” Exactly how does this *confirm* introgression?\n“Therefore, we used the statistic U20SAI, EAI, banteng (1%, 20%, and 100%) 20, which was\nequal to the number of SNPs within a genomic window where a particular allele was fixed at\na frequency of 100% in banteng but at a frequency less than 1% in SAI cattle or greater than\n20% in EAI cattle.” For both the banteng and gaur analyses you only have two individuals of\neach species. This means that you have essentially no ability to estimate allele frequencies\nin gaur and banteng. This means that all of your results from analyses of this type are\nsuspect, which in turn means that inference based on these results is suspect. There are\nmore than two of each of these publicly available in SRA and your analysis would be\nsignificantly strengthened by using what is available to increase your sample size, which will\nallow you to better estimate allele frequencies in these species.\nPaternal analysis\n“…(i) only present in at least two males but not in females; (ii) hemizygous…” Are these not\nrepresented as homozygous in the actual data? In reality, the males are hemizygous but\nthey are represented as homozygous when variant calling.\n“BEAGLE was used to impute missing alleles.” Were samples set to homozygous for\nimputation and you simply pretended as if this was an autosome? Exactly how was Beagle\nrun for Y-specific loci? The underlying model includes recombination rate, which is zero for\nthe Y. It’s unclear to me how this was performed.\nWhole mitogenome phylogeny\n“…samples with a depth of coverage lower than 100× were disregarded.” Why? The depth of\ncoverage on the MT is a function of the overall average genome coverage and the tissue\nsource for the DNA. This biases the analysis against samples that were sourced from semen.\nA more appropriate threshold would be to look at the average coverage of the MT relative to\nthe average of the autosomes. What you will find is that tissue sourced from semen will\nhave similar MT coverage to the autosomes while tissue sourced from anything other than\nsemen will have exceedingly high MT coverage. Regardless of coverage, you can assemble\nthe MT reads from all of the samples and then perform an evaluation of MT genome\ncompleteness versus coverage and set some (non-arbitrary) threshold for what samples you\nuse for downstream analysis. In summary, I think you may be leaving information on the\ntable by setting this arbitrary threshold.\nEstimation of effective population size and divergence time\nIt is unclear which variants were used and how the phasing was done. Both of these details\nneed to be fully documented since both impact the downstream analysis.\nSupplementary Fig. 4 & 5\nThis shows the results and the percentages for the first and second PCs are very small,\nespecially given the number of samples. Supplemental Fig 5A shows that the percentages\nget even smaller when PCA is performed on the indicine samples only. This data should be\nbetter filtered to get a more informative PCA. This raises further questions related to my\ncomments under the PCA section. I’m not convinced that this was done properly and\ntherefore the results and interpretation are suspect.\nMain manuscript which I’ll just refer to line numbers.\nL102 Change disposal to dispersal.\nLine 182 Fig.2. Personally, I hate these figures. They look pretty but show nothing of\nsubstantive value because of the resolution (genome coordinates). It is up to the authors as\nto whether to include it but I feel that it is a waste of space.\nL182-183 It is unclear why this is important or even relevant. Differentially expressed genes\nimplies that a regulatory variant has been selected differently between the two groups of\nanimals. However, a coding region variant may have been selected differently between the\ntwo groups and not change expression levels. What if you would have selected a different\nset of tissues? My point is that this analysis does nothing to further define what may be the\nunderlying cause of the region to be differentially selected between the two groups\n(assuming the region truly is a selection target) and simply adds noise to the discussion. I\nwould argue that this same observation could be made for a large number of regions\nrandomly selected from the genome. However, without quantifying this, you are making a\nlot of assumptions to relate gene expression differences to putative selection regions.\nL185 It appears to me that this is actually two regions. In fact, you list it as two regions but\nyou describe it as a single region. This needs to be reconciled.\nL192-198 All of this is speculation. While I am not opposed to some level of speculation, I\nam opposed to adaptive storytelling. If you truly identified a functionally plausible gene, and\nyou believe your underlying sequence data, then you should be able to identify a plausible\ncandidate mutation (or mutations) to explain your data. If you are going to state this, then\nwhy not do the follow-up analyses to try to identify the actual mutation(s)? I believe a deep\ndive into this would be far more valuable than what was done.\nL242-243 The SEMA3F Val650Ala variant appears to be 22:50162746 which is at frequency\n0.033 in the 1000 Bulls Run9 data and is at frequency of 0.034 in my own UMAG1 data\nbased on ~5500 samples. There are a lot of other protein altering variants in this gene, this\njust happens to be one of them. It is unclear how you arrived at this particular variant. I\nthink the manuscript would be significantly improved if the authors tried to dissect some of\nthese and do a more thorough analysis rather than simply list gene names.\nL262-265 I previously listed my concerns with this analysis. Adding to this, for almost all\nother analyses you used sliding windows. However, here you state that you used non-\noverlapping windows, why?\nL267 Sup. Figs 18-21. For Supp. Fig 19 has banteng in one part but gaur in the other. This\nneeds fixed.\nFor supp Fig 18 & 19, this looks like a lot of random noise to me and I suspect it has\nsomething to do with the banteng/gaur N=2 issue. Were consecutive windows merged? The\nsize of these regions should be informative for the timing of the introgression but nothing\nwas mentioned about this. I suspect due to the previous issues I’ve already raised.\nFor supp Fig 20 & 21, how was a p-value calculated? I could not find this.\nL281 Exactly how does this provide validation?\nL292-294 What is the call rate for these 11 loci in the raw data? Can you rule out any effect\nimputation may have had in this region? This is the most compelling evidence shown thus\nfar. I would recommend a deeper dive in this region to try to further strengthen your\ninference.\nL318 Supp Fig 45 legend says 2 panels but there is only 1. Fig 46 shows Ne increasing in\nrecent times which is contrary to everything we know about cattle demography. Perhaps\nthis relates to my next point.\nL325-344 This paragraph leads off with MSMC but it is unclear to me how you can use the Y\nand MT in an MSMC analysis when the underlying model is based on recombination?\nL402-446 This is all speculation. Again, some is useful but this is 1.5 pages. I would\nrecommend that a deeper dive into any one or two of these would be more valuable to the\nreader that speculating on all of them.\nL439-440 Proposes that EAI cattle may have introgression from a Bos-like ghost species as\nan explanation for hemoglobin-related genes. Since hybridization and introgression from\nknown Bos species is difficult, and these large numbers of divergent sequences appear in\nonly this family of genes, is it possible that this ’ghost’ group is a lost population of EAI-like\nIndicine cattle, or that these mutations were specific to the extant EAI clade without any\nintrogression?\nMethods\nL491-493 were duplicate reads marked or removed? The 1000 Bulls spec is for them to be\nmarked. I just want to be sure that what is stated is what was actually done. Along those\nlines, the spec has indel realignment and BQSR but that is not stated in the manuscript.\nPlease accurately specify what was actually done.\nL494 “…depth (for all individuals) > 1/3× and < 3×…” I have no idea what this represents,\nplease clarify.\nL500 Please specify BEAGLE parameters that were uses. If defaults were used, then state\nthat. Additionally, you specify SNPs here, were only SNPs used or did you also use indels.\nPlease specify. Did you make any attempt to evaluate the accuracy of imputation? If so, you\nshould state that. If not, I would encourage you to evaluate this and include this information\nin any filtering that you perform.\nL502 This is based on an annotation version, in which case you should specify the exact file\nor annotation version that was used to make it reproducible.\nL547 I’ve already discussed the multiple testing issue and choice of p<0.005.\nL557 Introgression analyses… How might a MAF threshold of 0.01 filtering affected these\nanalyses?\nL573 Same comment about the GTF version.\nL583 Already commented on the Y chromosome imputation.\nFigure 3 panels B-C and D-E appear to be switched relative to the legend. 10kb sliding\nwindow appears to be different than what is described in the M&M.\nI waive anonymity for all manuscripts and grants that I review and sign my reviews, Robert\nSchnabel."}, "author_response": "1 Response to Reviewers’ comments\n2 [Comment of Reviewer #1]\n5 genomes and 141 taurine cattle genomes from 57 and 17 populations, respectively. The quantity and\n6 quality of the material studied here, the methods used and the results obtained are of a high\n7 standard and deserve publication in a high-ranking journal. However, the main problem\n8 concerns the core zebu breed group and the outgroups. Regarding the core Zebu group: The\n10 are Gir, Nellore and Brahman, which were bred outside India for centuries in new environments\n11 and with sporadic interchange with Bos taurus breeds. Introgression of Bos taurus into the WGS\n12 of Gir, Nellore and Brahman used here is neither ruled out nor investigated by the authors.\n14 such as Sri Lanka, Pakisitan and Nepal (SAI group) use Bos taurus breeds to improve local\n15 cattle breeds. Therefore, the purity of any reference used for introgression studies should be\n16 tested first.\n18 description of the reference populations used for introgression studies.\n20 indicine (AMI) and were not included in the South Asian indicine (SAI) group. Cattle in Sir\n23 and pure indicine cattle that are included and verified in the results of admixture analysis. We\n24 provided this information in Supplementary Table 1. We selected indicine individuals without\n25 taurine ancestry. We apologize for not describing this in previous version of our manuscript.\n27 cattle have a pure genomic background. Our results showed that the reference SAI cattle did\n28 not carry either taurine or banteng/gaur ancestries.\n29\n30 [Comment of Reviewer #1]\n32 important introgression analyses in this study. Again, there is a possibility of introgression of\n33 cattle into Banteng and/or Baurs. This possibility has not been ruled out and has not been\n34 investigated. There are some unpublished results confirming cattle DNA segments in Bantengs.\n36 analysing cattle introgression in Yak, RFMix can produce false positive signals if the reference\n37 population is itself introgressed.\n40 been documented (Wu et al., 2018. Ref.15; Gao et al., DOI:10.1038/s41598-017-16438-7). To\n41 address your concern, we tested for introgression of cattle into banteng and gaurs. The f3\n42 statistic showed that there was no introgression into the banteng or gaur samples included in\n43 this study.\n45 and phylogenic trees of specific genes to verify this point, which in our view provided\n46 compelling evidence for introgression (Fig. 5D, Suppl. Figs. 16 to 28) and clustering of East\n48\n49 [Comment of Reviewer #1]\n51 analysis. I expect that the relatively small sample size of the outgroups (2 samples each) will\n52 affect the performance of RFMix and the determination of ancestral alleles and their frequencies\n53 in the U20 and U50 analyses. In fact, the original paper of RFMix already describes how the\n54 sample size of the reference panel affects its performance (https://www.ncbi.nlm.nih.gov/p\n55 mc/articles /PMC3738819 /figure/fig6/). How did the authors determine the robustness of their\n56 results?\n59 because of CITES regulations. However, we have added six banteng and three gaur samples to\n60 increase the number of samples of our reference group. And we reanalyzed the introgression\n61 from banteng or gaurs using U20 and RFMix analyse.\n63 gaur fragments and their homologous sequences in other bovine species belonging. This\n64 indicated unambiguously a banteng or gaur origin of the introgressed segments.\n65\n66 [Comment of Reviewer #1]\n68 well as incomplete lineage sorting. I recommend that the authors perform other independent\n69 approaches, e.g. those based solely on branch length and tree discordance, to validate these\n70 results.\n73 introgressed fragments that were likely caused by incomplete lineage sorting (ILS) according\n74 to a probability calculation (see Supplementary Note 4). We also analyzed the distribution of\n75 tree topologies of banteng or gaur fragments and their homologous sequences in other bovine\n76 species. This confirmed the banteng and gaur origins of the introgressions (Supplementary Fig.\n77 6).\n78\n79 [Comment of Reviewer #1]\n81 they should also include the matrix of residual values to show how integrating migration edges\n82 improved the % of variance explained.\n85\n86 [Comment of Reviewer #1]\n88 make sure that it is significantly different from the expected length because of its shared\n89 ancestral origin? Does the length of introgressed haplotypes tell us anything about the age of\n90 introgression?\n93 their introgressions. In our revised version, we modeled the expected length of ancestral\n94 sequence shared by indicine cattle, banteng and gaur on the basis of incomplete lineage sorting\n95 (ILS). The expected length of a shared ancestral sequence is L = 1/ (r × t), where r is the\n96 recombination rate of 1.23 cM/Mb (Weng et al. 2014), and t is the branch length between cattle\n97 and Asian wild bovine species (banteng and gaur). We used a generation time of 6 years and a\n98 divergence time of 1000 kya (Wu et al. 2018). We calculated the probability of ILS. The\n99 probability of a shared haplotype length derived from ILS according to the algorithm is 1 –\n102 probability < 0.05. We realize that this calculation, although according to the state-of-the-art,\n103 assumes a fixed recombination rate and does not account for the plausible effects of divergence\n104 between maternal and paternal haplotypes on the probability of recombinations. For this\n105 reason, we did not try to infer an age of the introgressions.\n106\n107 [Comment of Reviewer #1]\n109 stand-alone object, i.e. the readers should be able to understand it without reading other objects.\n111 message of this stand-alone object. This is not the case. I make a specific comment for one\n112 figure only (mentioned in the comments below), but it applies to all and is therefore a general\n113 comment.\n116\n117 [Comment of Reviewer #1]\n120 mentioned in the reference paper cited here.\n122\n123 [Comment of Reviewer #1]\n125 also into almost all south-eastern European cattle breeds up to the southern Alps\n126 (https://doi.org/10.1186/s12711-020-00560-8).\n129 well documented the male-mediated indicine admixture into African taurine cattle in the eastern,\n130 western and southern areas of the continent7,11and small but significant indicine introgression\n131 into almost all southeastern European cattle breeds12”.\n132\n133 [Comment of Reviewer #1]\n135 hypothesis about possible wild introgression in southern Indian zebu, the same is not mentioned\n136 for Gujarat, consider rewriting this sentence.\n138\n139 [Comment of Reviewer #1]\n141 the number of WGS used for each species separately.\n143 sequences used for each species as follows in the revised version: “We also used sequencing\n144 data of 22 whole genomes from six other bovine species, including two bison, two wisent, five\n145 gaur, eight banteng, two yak, and two water buffaloes, as outgroups or for introgression\n146 analysis.\n147 [Comment of Reviewer #1]\n152 genetically in an intermediate position between SEAI and EAI”. According to the comments of\n154 ancestries. The PCA and phylogenetic tree almost completely separated the three indicine\n155 geographic lineages of SAI, AFI, and EAI. SWCI was in an intermediate position between SEAI\n156 and EAI.\n157\n158 [Comment of Reviewer #1]\n160 combined the three indicine groups (SAI, AFI, EAI) and left the other two (SEAI and SWCI)\n161 for this analysis?\n164 two groups. Thus, we believe that the three indicine groups of SAI, AFI, and EAI adequately\n165 represent the indicine cattle ancestry.\n166\n167 [Comment of Reviewer #1]\n169 groups of animals/populations: 1) taurine and 2) SAI, AFI and EAI indicine populations. Both\n170 groups lived in environments that differed by a variety of factors, including temperature. The\n171 authors found numerous significantly differentiated regions and attribute them all to adaptation\n172 to high temperatures. The idea of reverse genomics is to find the outliers and then look for\n173 possible causes. Only after some causality analysis can we conclude that some of the outliers\n174 are an adaptation to high temperatures or something else. The authors should consider\n175 significant selection signatures as anonymous candidates. These anonymous candidates are\n176 subject to positive selection into a group of animals that share a specific environment. These\n177 anonymous candidates could, but need not, be caused by one specific environmental factor.\n180 temperature. We have modified the headings of Lines 174 and 175 as follows:\n181 “The ancestral adaptation of indicine cattle\n184 warranted to test their role in heat adaptation or other differences between indicine and taurine\n185 cattle.”\n186\n187 [Comment of Reviewer #1]\n189 gender and age) are quite heterogenous, I will take this result with pinch of salt, moreover, how\n190 many other genes were differentially expressed and not under positive selection?\n193\n194 [Comment of Reviewer #1]\n196 indicine cattle by controlling hair length and/or thickness. However, we (readers) are not\n197 informed whether indicine and taurine cattle differ in hair length and thickness? I cannot\n198 find any such information in this manuscript.\n200 between taurine and indicine cattle in the Introduction as follows:\n201 “Indicine cattle are recognized by their thoracic hump, low metabolic rate, many large sweat\n202 glands, large skin surface, and short smooth coat 5.”\n203\n204 [Comment of Reviewer #1]\n206 broad area), but all outliers are defined as candidates for adaptation to semi-arid environmental\n207 stress.\n209 not necessarily indicate a specific adaptation. We have modified the headings of Lines 174 and\n210 175:\n211 “The ancestral adaptation of indicine cattle\n214 their roles in heat adaptation or other differences between indicine and taurine cattle”.\n215\n216 [Comment of Reviewer #1]\n218 therefore, I would remove it or declare as confirmation.\n220\n221 [Comment of Reviewer #1]\n224 emerged as a new haplogroup during the indicine westward migration, but as a haplogroup that\n225 was already present in the area before the indicine introgression. The authors should consider\n226 both the time of coalescence and the time of migration when drawing their conclusions.\n229 taurine cattle that were already present in Africa prior to the introduction of indicine cattle (see\n230 lines 455-458; ref. 7 (Kim et al., 2020, The mosaic genome of indigenous African cattle as a\n231 unique genetic resource for African 691 pastoralism. Nature Genetics 52, 1099-1110); and ref.\n232 11 (Kim et al., 2017, The genome landscape of indigenous African cattle. Genome Biology 18,\n233 34).\n234\n235 [Comment of Reviewer #1]\n237 could be a reason why the west-east cline is missing and an abrupt transition is observed. I did\n238 not understand this section of the text. There are several problems in the text itself and in the\n239 corresponding figures. The legends of the figures are generally not informative enough. Let's\n240 just take Supplementary Figure 17 as an example. The reader cannot understand this figure from\n241 its legend. Even more, there are many trees labelled with the same A, B, C,.... and we do not\n242 know which tree represents which haplotype, etc. I also do not know what to conclude from\n243 figure 50, i.e. the legend does not help me to understand the message from line 355-357.\n246 figures to ensure that they are sufficiently informative.\n248 species belonging to the bovine species to confirm banteng and gaur introgression.\n250 provides a geographic contour map of banteng/gaur introgression proportions in 16 East Asian\n251 indicine (EAI) breeds. The proportions of banteng and gaur introgressions were calculated by\n253 visualized using the ArcMap component of the ArcGIS software suite. EAI cattle in the\n254 southeastern coastal region of China show the highest level of banteng and gaur ancestries. We\n255 also modified the legends of Supplemental figure 50 (now is Supplementary Fig. 29) in order to\n256 explain how the geography with three large rivers in a mountainous regions impede the gene\n257 flow between SEAI and SWCI.\n258\n259 [Comment of Reviewer #1]\n262\n263 [Comment of Reviewer #1]\n265 cattle could be helpful here (DOI: 10.1111/eva.13315 and doi: 10.1111/age.13104).\n267 mitochondrial and Y-chromosomal variations in cattle. We now refer the first paper on taurine\n268 mtDNA in line 124. Since the second paper of Escouflaire et al. (2021) described Y-\n269 chromosomal variation in French taurine cattle, we do not refer it due to its irrelevance to our\n270 analyses.\n271\n272 [Comment of Reviewer #1]\n275 the spread of domestic sheep in Europe. However, for sheep, there is archaeological evidence\n276 for a sheep-less Europe before the Neolithic. Do you have any evidence that East Asia was a\n277 cattle-less region before the spread of Bos indicus? Africa was not cattle-less either, and as a\n278 result of paternal introgression we only observe haplogroup T. Please clarify.\n280 cattle prior to the arrival of indicine DNA (Felius et al., 2014, doi:10.3390/d6040705; Zhang\n281 et al., 2018; doi.org/10.1186/s12863-018-0705-9). Approximately 20% of South Chinese cattle\n282 still contain taurine mtDNA with high diversity (Gao et al., DOI:10.1038/s41598-017-16438-\n283 7; Xia et al., 2018; doi: 10.1111/age.12749). This is now indicated in the text as follows:\n284 “Indicine cattle may have entered East Asia between 3,500 and 2,500 YBP well after the arrival\n285 of taurine cattle (ref. 10, Felius et al., 2014, doi:10.3390/d6040705; Zhang et al., 2018;\n286 doi.org/10.1186/s12863-018-0705-9; Gao et al., doi:10.1038/s41598-017-16438-7; Xia et al.,\n287 2018; doi: 10.1111/age.12749)”.\n288\n289 [Comment of Reviewer #1]\n291 trading routes from South to East Asia. Therefore, the integrated genetic analyses as mentioned\n292 in the paper may have difficulty supporting inland trading even if this is correct.\n295 autosomal data support a coastal route for the first migration wave to Southeast Asia as the\n296 main entry point of indicine cattle into East Asia”.\n298 as the rivers and mountainous areas in Southwest China, which does not favor an inland\n299 migration of cattle (Supplemental Figure 39).\n300\n301 [Comment of Reviewer #1]\n305\n306 [Comment of Reviewer #1]\n309 “SriLanka” exists.\n312 name of Sri Lankan cattle.\n313\n314 [Comment of Reviewer #1]\n317 2018 and should be from India. We know that India is very restrictive when it comes to foreign\n318 use of biological resources. There is no Indian collaborator in Chen et al. 2018, so I do not see\n319 any legal way to analyze Indian biological material outside India or without an Indian\n320 collaborator. Was this Hariana cattle sampled in India or elsewhere? Please clarify and correct\n321 if necessary in Supplementary Table 2.\n324 animals have been in his collection since the early 1990s, and they were sampled at the Indian\n326 first paper including these samples (Anim Genet. 1994; 25(4): 265-71. doi: 10.1111/j.1365-\n327 2052.1994.tb00203.x.) did have a collaborator, D S Balain, from India.\n328\n330 [Comment of Reviewer #2]\n332 not been previously characterized. Overall, I find the manuscript well written, although it\n333 appears to have been done in parts by different groups of individuals with different writing\n334 styles and slightly different nomenclature. This is inevitable with a project of this size. I would\n335 encourage the corresponding authors to go through the entire manuscript to make it more\n336 uniform. I do have a number of concerns related to the details that were presented, or not\n337 presented, that make a final determination not possible at this time. I always start my reviews\n338 with the materials and methods but since there was such a substantial amount of information in\n339 the supplemental information, I will start there. Unfortunately, the supplemental document does\n340 not have line nor page numbers, so I’ll quote specific lines.\n341\n342 [Comment of Reviewer #2]\n346 because these are run IDs. You should use the BioSample ID. I downloaded this table and\n347 compared these SRR to the NCBI database and you have two duplicate samples. SRR2016752\n348 and SRR2016754 are the same SAMN3387026 individual. With only a single duplicate, it is\n349 unlikely that any of the results would have been significantly affected. However, this table needs\n350 to be corrected to present the BioSample and if data reanalysis is performed, one of these two\n351 should be removed.\n353 the results. We also provided the BioSample ID for all samples in Supplementary Table 1 and\n354 removed one duplicated sample.\n355\n356 [Comment of Reviewer #2]\n358 “Runs of homozygosity (ROH) were identified using the --homozyg option implemented in\n360 genome.” This is a very small window considering the density of markers. It is well recognized\n361 that results are sensitive to settings for this analysis (see https://doi.org/10.1186/s12864-020-\n362 6463-x). I question whether or not the settings used were appropriate. Specifically, by using too\n363 small of a window size and other inappropriate settings, far too many ROH are detected.\n365 samples that have estimated ROH totaling greater than 1/3 of the genome and a large number\n366 of samples with ROH greater than half the genome. This is unrealistic. This analysis needs to\n367 be redone with appropriate settings, which may need to be determined empirically given the\n368 data. It would also be useful at this time to examine any correlation\n369 between average genome coverage and number/length of ROH as one would predict that lower\n370 coverage samples will have heterozygous sites undercalled which may manifest as ROH when\n371 in fact there are not ROH if heterozygous positions were called accurately.\n374 coverage < 90% and used 331 individuals for ROH analysis. We then used imputed SNPs to\n375 detect ROHs with a minimum length of 100 kb and containing at least 50 SNPs using PLINK\n376 v1.9 software. Additionally, three heterozygous SNPs were allowed per ROH.\n378 varying this setting from 10 to 200 SNPs (step size of 10 SNPs). The other parameters were set\n379 to a minimum density threshold (50 SNPs), a large gap (1000 kb), a minimum length (50 kb), a\n380 maximum number heterozygous SNPs in a scanning window (3), and a scanning window\n381 threshold level (0.05). The results suggested that increasing scanning window size led to the\n382 decrease in both number and total length of the estimated ROHs (Fig. 1), however, there are\n383 still too many ROH; When the scanning window size was 100 and 150 Mb, the largest numbers\n384 of ROHs were 7914 and 6735, respectively.\n385\n386\n388 for all individuals from each cattle geographic group. Each dot represents an individual.\n390 to filter the small ROHs. The other parameters were set to a minimum density threshold (50\n392 and a scanning window threshold level (0.05). The results show that our settings get the\n393 expected number (maximum number is 4570) and total length (maximum length is 1,423,160\n395 these results in the Supplementary Fig. 2. And these results shows that the level of inbreeding\n396 measured by ROHs was lower in indicine cattle than in taurine cattle.\n397\n399 scanning window size is 50 SNPs. (B) the scanning window size is 100 SNPs. The minimum ROH length\n400 was set to 100 Kb for excluding short ROHs.\n402 validity of the ROH analysis (Meyermans et al. 2020). The scanning window size setting (--\n403 homozyg-window-snp) was investigated by varying this setting from 10 to 200 SNPs (step size\n404 of 10 SNPs). The unchanged parameter set to a minimum density threshold (50 SNPs), a large\n405 gap (1000 kb), a minimum length (50 kb), a maximum number heterozygous SNPs in scanning\n406 window (3) and a scanning window threshold level (0.05). Consequently, genome coverage was\n407 higher than 99% for all breeds, which means that the given settings allowed ROH detection for\n408 more than 99% of all autosomes (Fig. 3).\n409\nROH\n411 breed/populations. (A to C) three taurine cattle breeds (Jersey, Mongolian, and Yanbian). (D to F) three\n412 indicine cattle breeds (Longlin, Weizhou, and Red Sindhi).\n413\n414 [Comment of Reviewer #2:]\n416 “We removed \"all LD\" using…” First, please rephrase this. You did not remove “all LD”. I have\n417 the same objection here to the parameters used. First, this is a very small window and probably\n418 left a significant number of loci with r^2 > 0.20 which are more than 50kb apart. Given\n419 everything that we know regarding LD in bovine genomes, you should empirically find the\n420 appropriate setting for this analysis to achieve the desired thinning. This does not have to be\n421 done genome-wide and one could simply use 10Mb from any autosome to estimate what the\n422 appropriate parameters should be. My recommendation would be to use chr25 since we know\n423 there are no assembly issues. For a 10Mb region, calculate *all* the pairwise LD values using\n424 ‘-r2 yes really’ and examine the distribution of these values. This figure can be included in the\n425 supplement to provide justification for the settings chosen. As the manuscript is currently\n426 written, the settings are arbitrary and likely impacted the results and interpretation.\n429 provided these results in Supplementary Fig. 4. The half value of LD decay (r2) is 0.11, so we\n430 performed LD-based pruning for the genotype data using PLINK v1.9 with ‘-indep-pairwise 50\n431 10 0.1’.\n432\n434 half value of LD decay is 0.13.\n435\n437 20 10 0.1’ options to perform LD based pruning for the genotype data and used these data for\n439 admixture analysis (Supplementary Fig. 5).\n441 indep-pairwise 50 10 0.1’ option for LD pruning.\n442\n443\n445 against PC2 (A) and PC1 against PC3 (B). A total of 2,996,368 LD pruned SNPs were used for\n447\n448\n450 against PC2 and PC1 against PC3. A total of 7,192,063 LD pruned SNPs were used for PCA\n451 with the parameter ‘-indep-pairwise 20 10 0.1’.\n452\n453 [Comment of Reviewer #2:]\n454 “PCA was performed on the genome-wide unlinked SNP dataset…” Please change the word\n455 “unlinked” throughout the manuscript. You do not have any unlinked variants. Perhaps use “LD\n456 pruned” as a replacement.\n458 throughout this revised version of our manuscript.\n459\n460 [Comment of Reviewer #2]\n462 “…200 bootstrap replicates were performed…” Why 200? Why not 10, or 10,000. Is the default\n463 of 200 appropriate for your data? “…used to determine the optimal ancestry number…” Exactly\n464 how did you determine the optimal number? There is quite a bit of literature on this, but nothing\n465 was stated or cited.\n468 default setting, so bootstrap replicates were not used. We used the delta K method to choose the\n469 optimal K, and we added the K values to the Supplementary Table 4. And a sentence is added\n470 as follows:\n471 “ADMIXTURE v.1.3.0 was used to quantify genome-wide admixture among cattle\n472 breeds/populations 82 and run for each possible group number (K = 2 to 8), where K is the\n473 assumed number of ancestral populations. The delta K method was used to choose the optimal\n475 [Comment of Reviewer #2]\n477 “…A window size of 1,000 SNPs was used to account for linkage disequilibrium…” This\n478 should be unnecessary if you properly LD pruned in the first place. This brings up another issue\n480 which set of loci you are using for which analysis. In this section, it appears that you are using\n481 the “raw” 67M variants. However, it seems to me that if the objective is to “identify closely\n482 related individuals” then the LD pruned dataset from the PCA section would be appropriate.\n484 be to provide a supplemental table that lists the major analyses performed and exactly which\n485 data manipulations were done for each analysis. Since data manipulations may impact one\n486 analysis but not another, it would be helpful to the reader to understand exactly what was done\n487 and which analyses used the same loci.\n491 and ML trees as follows: “To identify relationships among individual cattle, a total of\n492 67,162,108 autosomal SNPs were used to construct a NJ tree with PLINK v.1.9 based on the\n493 matrix of pairwise genetic distances 6 (Fig. 1). FigTree v.1.4.3 10 was used to visualize the NJ\n494 tree. Then, we inferred a population-level phylogeny using the ML approach implemented in\n496 using the --indep-pairwise 50 5 0.1 option of PLINK v.1.9 6. A total of 15,228,801 LD pruned\n500\n501 [Comment of Reviewer #2]\n503 “…were calculated for 50-kb windows with a 20 kb step across the autosomes…” Line 907 of\n504 the main manuscript says that a 100kb window was used, which is it? I have issues with what\n505 is specified here. Typically, these types of analyses are performed where the window size is an\n506 integer multiple of the step size. For instance 25 kb step and 50 kb window or 20kb step and\n507 100kb window. This needs to be addressed throughout the manuscript and I would recommend\n508 using the same parameters for all analyses that use this type of setting (where appropriate). This\n509 will make comparisons between analyses much easier.\nST\n512 sliding window approach with windows of 50 kb and a step size of 20 kb.\n514 kb and a step size of 20 kb.\nST\n516 taurine), F , XP-EHH, and PBS, we used 50 kb and a step size of 20 kb. For U20 and U50\nST\n517 analyses, we reanalyzed the data and used a 50 kb window and a step size of 20 kb too.\n518\n519 [Comment of Reviewer #2]\n520 “…these windows harbored 117 candidate genes in indicine cattle…” I don’t think these 117\n521 genes are specific to indicine cattle and you can probably just say cattle. By saying indicine\n522 cattle here one might think you are implying that these are not in taurine cattle, which isn’t\n523 the case.\n525 the reference to indicine cattle. The sentence is revised as follows: “These windows harbored\n526 117 candidate genes.”\n527\n528 [Comment of Reviewer #2]\n529 “…the expression levels of these 117 candidate genes were validated in nine different tissues\n530 of taurine and indicine cattle.” I’m not sure exactly what you are validating? The fact that there\n531 are RNA-seq data in those tissues for these annotated genes simply validates that the annotation\n532 is correct.\n534 this part from the manuscript.\n535\n536 [Comment of Reviewer #2]\n537 “Candidate genes under selection were defined as those overlapped by sweep regions or within\n538 20 kb of the signals.” 20kb seems arbitrary. What justification do you have for choosing 20kb?\n540 any overlap. It seems to me that allowing for a gene to be outside your sweep region, by some\n541 arbitrary distance, speaks to your confidence in defining the boundaries of sweep regions. Some\n542 justification should be given for whatever value you choose.\n545 taken from the selected regions identified with the three methods, without extending the signal\n546 by 20 kb. This sentence has been revised as follows: “The candidate genes selected in all\n547 indicine cattle were defined as the genes with overlapped signals in any two of the three\n548 selection methods (θπ ratio (indicine/taurine), F , and XP-EHH)”.\nST\n549\n550 [Comment of Reviewer #2]\n551 “…were plotted using a 10-kb sliding window…” Again, arbitrary. Why not use the same\n552 window size for plotting that you used for detection?\n556 use the detection window (50 kb), there are too few SNPs to draw the graph, and the real signal\n557 for the target region will be overlooked.\n558\n559 [Comment of Reviewer #2]\n561 of 20 kb…” Same comment on windows.\n564 window with a 20 kb step size across the autosomes, including CLR, F , and XP-EHH.\nST\n565\n566 [Comment of Reviewer #2]\n567 “…the overlap windows of P values less than 0.005 (Z test) of each method were considered\n568 candidate signatures of selection…” Why 0.005? Again, some justification is warranted. It\n569 seems to me that you have an extraordinarily large number of windows that were tested which\n570 implies that some adjustment needs to be considered for multiple testing. Appropriately\n571 accounting for multiple testing here could significantly change the detected regions which in\n572 turn could significantly impact all downstream inference based on these results. As written, it\n573 is impossible to determine exactly what was done and what impact it may have on the results\n574 or interpretation.\n577 and were further corrected by multiple testing using the Benjamin-Hochberg false discovery\n578 rate (FDR) method. We hope this justified.\n579\n580 [Comment of Reviewer #2]\n581 “Only pathways or annotations with a Bonferroni-corrected P < 10−2 were retained\n582 (Supplementary Tables 7-10).” This should probably be written as P<0.01 to be consistent\n583 throughout the manuscript.\n585 throughout the manuscript.\n586\n587 [Comment of Reviewer #2]\n589 “Phylogenetic analyses of these segments confirmed the banteng or gaur introgression into\n590 specific EAI individuals.” Exactly how does this *confirm* introgression? “Therefore, we used\n591 the statistic U20SAI, EAI, banteng (1%, 20%, and 100%) 20, which was equal to the number\n592 of SNPs within a genomic window where a particular allele was fixed at a frequency of 100%\n593 in banteng but at a frequency less than 1% in SAI cattle or greater than 20% in EAI cattle.” For\n594 both the banteng and gaur analyses you only have two individuals of each species. This means\n595 that you have essentially no ability to estimate allele frequencies in gaur and banteng. This\n596 means that all of your results from analyses of this type are suspect, which in turn means that\n597 inference based on these results is suspect. There are more than two of each of these publicly\n598 available in SRA and your analysis would be significantly strengthened by using what is\n599 available to increase your sample size, which will allow you to better estimate allele frequencies\n600 in these species.\n603 to 27) and support that introgressed segments of East Asian haplotypes were clustered within\n604 bangteng or gaur haplotype groups. According to your and reviewer 1’s suggestions, we have\n605 added six banteng samples and three gaur samples, and we recalculated the U20 and U50\n606 statistics. We have added the number of whole genome sequences used for each species as\n607 follows in the revised version: “We also used sequencing data of 22 whole genomes from six\n608 other bovine species, including two bison, two wisent, five gaur, eight banteng, two yak, and\n609 two water buffaloes, as outgroups or for introgression analysis.”\n610\n611 [Comment of Reviewer #2]\n614 hemizygous but they are represented as homozygous when variant calling.\n617 follows: “Only the SNPs called in the MSY region that met the following criteria were retained:\n618 (1) present in at least two males but not in females and (2) no heterozygous site.”\n619\n620 [Comment of Reviewer #2]\n621 “BEAGLE was used to impute missing alleles.” Were samples set to homozygous for\n622 imputation and you simply pretended as if this was an autosome? Exactly how was Beagle run\n623 for Y-specific loci? The underlying model includes recombination rate, which is zero for the Y.\n626 assumption that they may be belonged to autosomes. In this study, we used imputed data only\n627 to construct a median-joining (MJ) network. To avoid misleading results, we reconstructed the\n628 network using unimputed data. We deleted this imputation-related part to avoid misleading.\n629\n630 [Comment of Reviewer #2]\n632 “…samples with a depth of coverage lower than 100× were disregarded.” Why? The depth of\n633 coverage on the MT is a function of the overall average genome coverage and the tissue source\n634 for the DNA. This biases the analysis against samples that were sourced from semen. A more\n635 appropriate threshold would be to look at the average coverage of the MT relative to the average\n636 of the autosomes. What you will find is that tissue sourced from semen will have similar MT\n637 coverage to the autosomes while tissue sourced from anything other than semen will have\n638 exceedingly high MT coverage. Regardless of coverage, you can assemble the MT reads from\n639 all of the samples and then perform an evaluation of MT genome completeness versus coverage\n640 and set some (non-arbitrary) threshold for what samples you use for downstream analysis. In\n641 summary, I think you may be leaving information on the table by setting this arbitrary threshold.\n644 semen will have MT coverage similar to that of the autosomes, which will not be sufficient to\n645 assemble full mitogenomes. We apologize for our unclear description. We used stricter criteria\n646 for mitogenomes, and we have revised and supplemented our methods.\n648 first selected all indicine cattle in our dataset for mitogenome analysis, and then we selected\n649 only mitogenomes that were successfully assembled by MIA software and filtered mitogenomes\n650 with a gap length > 1 bp. According to these criteria, we compared our results with those of a\n651 previous study. We used more criteria, and an additional 13 samples were filtered. A total of\n652 329 mitogenomes assembled in this study and 18 reference mitogenomes were used for the final\n653 analysis. We reanalyzed all mitogenomes in the manuscript. The results were similar to our\n654 previous results.\n655\n656 [Comment of Reviewer #2]\n659 need to be fully documented since both impact the downstream analysis.\n662 time.\n664 (MSMC2) method was used to model the population history of the three core indicine groups\n665 (EAI, SAI, and AFI) and to infer historical changes in their effective population size and\n666 population separation. We applied this method to all groups with two deep-coverage (>14 ×)\n667 individuals per group. All sample sets of filtered variant calls were used for imputation and\n668 phasing using Beagle v4.1 with default parameters 4, except for “niterations” which was set to\n669 10. For each individual, DR2 value in INFO column of the “phase.vcf” file was used to filter\n671 recommended in the documentation of the software. For the calculation of effective population\n672 size, the parameter of MSMC2 was set to “msmc2 -t 10 -p 1*2+25*1+1*2 -I 0,1,2,3” and\n673 “msmc2 -t 10 -p 1*2+25*1+1*2 -I 4,5,6,7”. For the calculation of population separation, the\n674 parameter of MSMC2 was set to “msmc2 -t 8 -P 0,0,0,0,1,1,1,1 -s -p 1*2+25*1+1*2”. For\n675 effective population size inference, two individuals (4 phased haplotypes) from each population\n676 were used. The time scale in generation time at g = 6 and a mutation rate per generation at μg\n677 = 1.26×10-8 were used.”\n678\n679 [Comment of Reviewer #2]\n682 given the number of samples. Supplemental Fig 5A shows that the percentages get even smaller\n683 when PCA is performed on the indicine samples only. This data should be better filtered to get\n684 a more informative PCA. This raises further questions related to my comments under the PCA\n685 section. I’m not convinced that this was done properly and therefore the results and\n686 interpretation are suspect.\n688 performed PCA again using new LD pruned data with the parameters “-indep-pairwise 50 10\n689 0.1”. The new results are shown in Supplementary Figs. 4 and 5.\n690\n691 [Comment of Reviewer #2]\n695\n696 [Comment of Reviewer #2]\n698 substantive value because of the resolution (genome coordinates). It is up to the authors as to\n699 whether to include it but I feel that it is a waste of space.\n701 the Supplementary Fig. 9.\n702\n703 [Comment of Reviewer #2]\n705 implies that a regulatory variant has been selected differently between the two groups of animals.\n707 and not change expression levels. What if you would have selected a different set of tissues?\n709 of the region to be differentially selected between the two groups (assuming the region truly is\n710 a selection target) and simply adds noise to the discussion. I would argue that this same\n711 observation could be made for a large number of regions randomly selected from the genome.\n713 expression differences to putative selection regions.\n715\n716 [Comment of Reviewer #2]\n718 describe it as a single region. This needs to be reconciled.\n722\n723 [Comment of Reviewer #2]\n725 opposed to adaptive storytelling. If you truly identified a functionally plausible gene, and\n726 you believe your underlying sequence data, then you should be able to identify a plausible\n727 candidate mutation (or mutations) to explain your data. If you are going to state this, then\n728 why not do the follow-up analyses to try to identify the actual mutation(s)? I believe a deep\n729 dive into this would be far more valuable than what was done.\n731 tentative, and we have adapted our text as follows:\n732 “Seven of the other 75 genes located in the topmost significant sweeps are functionally\n733 associated with heart development, blood circulation, DNA damage, and light response.\n735 indicine and taurine cattle.”\n736\n737 [Comment of Reviewer #2]\n739 0.033 in the 1000 Bulls Run9 data and is at frequency of 0.034 in my own UMAG1 data based\n740 on ~5500 samples. There are a lot of other protein altering variants in this gene, this just happens\n741 to be one of them. It is unclear how you arrived at this particular variant. I think the manuscript\n742 would be significantly improved if the authors tried to dissect some of these and do a more\n743 thorough analysis rather than simply list gene names.\n745\n746 [Comment of Reviewer #2]\n748 other analyses you used sliding windows. However, here you state that you used non-overlapping\n749 windows, why?\n751 of other methods, we used the same standards (50 kb window and 20 kb step size) for\n752 recalculation in the revised manuscript.\n753\n754 [Comment of Reviewer #2]\n756 fixed. For supp Fig 18 & 19, this looks like a lot of random noise to me and I suspect it has\n757 something to do with the banteng/gaur N=2 issue. Were consecutive windows merged? The\n758 size of these regions should be informative for the timing of the introgression but nothing was\n759 mentioned about this. I suspect due to the previous issues I’ve already raised.\n763 five gaur samples, and we merged consecutive windows. As we used only the results of the U20\n764 statistic, we plotted only adaptive introgressed regions with higher frequencies; please see the\n765 results.\n766\n767 [Comment of Reviewer #2]\n770 transformed values using the standard normal distribution and were further corrected by\n771 multiple testing using the Benjamin-Hochberg false discovery rate (FDR) method. We have\n772 added this information to the manuscript.\n773\n774 [Comment of Reviewer #2]\n779 indicine cattle, and other bovine species clustered East Asian indicine cattle with banteng or\n780 gaur, thus confirming banteng or gaur introgression into the East Asian indicine genomes.\n781\n782 [Comment of Reviewer #2]\n784 imputation may have had in this region? This is the most compelling evidence shown thus\n785 far. I would recommend a deeper dive in this region to try to further strengthen your\n786 inference.\n789 and the missing rate is 0.01 to 0.21.\nChromosome Position Alleles Number of missing alleles Missing rate\n25 216559 160 0 0\n25 216571 160 0 0\n25 216581 160 0 0\n25 216850 160 16 0.1\n25 219613 160 2 0.0125\n25 219625 160 4 0.025\n25 219635 160 10 0.0625\n25 222783 160 0 0\n25 222849 160 0 0\n790\n791 [Comment of Reviewer #2]\n793 recent times which is contrary to everything we know about cattle demography. Perhaps this\n794 relates to my next point.\n796 in Supplementary Fig. 45 (now Supplementary Fig. 33). For Supplementary Fig. 46 (now\n798 associated with postdomestication expansion, as observed in other domesticated species, and\n799 with a more recent diffusion of I1a within East Asia.\n800\n801 [Comment of Reviewer #2]\n803 and MT in an MSMC analysis when the underlying model is based on recombination?\n807 the paragraph to avoid misunderstanding: “Using an empirical Bayesian approach with\n809\n810 [Comment of Reviewer #2]\n812 that a deeper dive into any one or two of these would be more valuable to the reader that\n813 speculating on all of them.\n815 speculation. We agree that this text is too detailed for a discussion of our tentative results. It\n816 has been condensed to 17 lines. We also made the preceding part of the Discussion more to the\n817 point.\n818\n819 [Comment of Reviewer #2]\n821 explanation for hemoglobin-related genes. Since hybridization and introgression from known\n823 family of genes, is it possible that this ’ghost’ group is a lost population of EAI-like\n825 introgression?\n827 the ‘ghost’ species is not compatible with an EAI origin but instead indicates a sister species of\n828 gaur. We also used an ancient kouprey genome (2× coverage) and gayal genomes to genotype\n829 this region. The polygenetic tree showed that EAI was influenced not only by banteng and gaur\n830 but by gayal and extinct kouprey species as well. Considering the low coverage of kouprey and\n831 the hybrid ancestry of gayal, we did not analyze the gene flow between these two species and\n833\n834 [Comment of Reviewer #2]\n837 marked. I just want to be sure that what is stated is what was actually done. Along those lines,\n838 the spec has indel realignment and BQSR but that is not stated in the manuscript. Please\n839 accurately specify what was actually done.\n842 map reads and we did not use indel realignment. We have revised the method in the manuscript.\n844 (http://www.1000bullgenomes.com/) (Supplementary Note 1). We removed low-quality bases\n845 and artifact sequences using Trimommatic v.0.39, and all clean reads were mapped to the cattle\n846 reference assembly (ARS-UCD1.2) and Btau_5.0.1 Y BWA-MEM (v.0.7.13-r1126) with default\n847 parameters 1. We then used SAMtools v.1.9 to sort bam files. For the mapped reads, potential\n849 (http://broadinstitute.github.io/picard). ‘BaseRecalibrator’ and ‘PrintReads’ of the Genome\n851 recalibration (BQSR) with the known variant file (ARS1.2PlusY_BQSR_v3.vcf.gz) provided by\n852 the 1000 Bull Genomes Project.\n855 ‘SelectVariants’, respectively. To avoid possible false-positive calls, we used VariantFiltration\n856 as recommended: (1) SNP clusters with ‘-clusterSize 3’ and ‘-clusterWindowSize 10’ options;\n857 (2) SNPs with mean depth (for all samples) < 1/3× and > 3× (×, overall mean sequencing\n858 depth across all samples); (3) quality by depth, QD < 2; (4) phred-scaled variant quality score,\n861 sum test, ReadPosRankSum < −8 were filtered. We then filtered out nonbiallelic SNPs and SNPs\n862 with missing genotype rates > 0.1. The whole genome sequencing data from six other bovine\n863 species were mapped in the same way. We used the 67,162,108 SNPs as a reference list to\n864 genotype the combined set of 495 cattle samples and 22 samples of six other bovine species,\n865 resulting in 67,145,163 SNP data with wild species data. The two final SNPs genotyping\n866 datasets were imputed and phased using BEAGLE v.4.0 3 with default parameters and filtered\n867 by DR2 < 0.9 (Supplementary Table 3). The remaining SNPs were annotated according to their\n868 positions using SnpEff v.4.3 4. We also summarized samples and SNPs used for different\n869 analyses in Supplementary Table 3.\n870\n871 [Comment of Reviewer #2]\n873 please clarify.\n875 mean sequencing depth (over all samples) < 1/3× and > 3× (×, overall mean sequencing depth\n876 across all samples) were removed”.\n877\n878 [Comment of Reviewer #2]\n880 that. Additionally, you specify SNPs here, were only SNPs used or did you also use indels.\n882 should state that. If not, I would encourage you to evaluate this and include this information\n883 in any filtering that you perform.\n885 of BEAGLE v.4.0 for imputation. We also detected accuracy according to DR2. In this version,\n886 we filtered sites with DR2 < 0.9, resulting in a total of 65,336,403 SNPs (Supplementary Table\n887 3).\n888\n889 [Comment of Reviewer #2]\n891 or annotation version that was used to make it reproducible.\n893 assembly GCF_002263795.1 for autosome annotation and RefSeq assembly\n895\n896 [Comment of Reviewer #2]\n900 and were further corrected by multiple testing using the Benjamin-Hochberg false discovery\n901 rate (FDR) method.\n902\n903 [Comment of Reviewer #2]\n905 analyses?\n907 detection of rare alleles of wild bovine species.\n908\n909 [Comment of Reviewer #2]\n914\n915 [Comment of Reviewer #2]\n918 version of manuscript.\n919\n920 [Comment of Reviewer #2]\n922 window appears to be different than what is described in the M&M.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nGeneral comments\nThe main message of this manuscript is that East Asian Indicine (EAI) cattle are descended\nfrom pure Indicine cattle originating in the South Asian core area (i.e. SAI, Indus Valley).\nFollowing the introduction of these pure Indicine animals via coastal routes, rapid and\nsuccessful adaptation to hot and humid environments occurred through introgression of\nBanteng and/or Gaur genes. The age of introduction of the pure Indicine cattle to the new\nenvironment is estimated to be about 3000 years, i.e. 600 generations of cattle. That mean\nthe first contact to and the introgression of Banteng and/or Gaur should be about 600\ngenerations ago or younger. These results are supported by cumulative evidence from\ndifferent analyses.\nAdaptive introgression is a very important evolutionary mechanism, but currently it is\nbecoming more popular and many signals are interpreted as introgression, often without an\nalternative hypothesis.\nMy main comments in the first report were about the quality of the reference groups, which\nare important for the above mentioned main message of this manuscript. I asked questions\nlike:\n1) Do the authors ensure that the SAIs are pure indicines without a genomic taurine\nsegment?\n2) Do the authors ensure that the Banteng samples are pure and do not contain a genomic\nsegment from indicine or taurine cattle?\n3) Do the authors ensure that the Gaur samples are pure and do not contain a genomic\nsegment from indicine or taurine cattle?\nIn his rebuttal letter, the authors state that this is ensured, but there are some\ninconsistencies that worry me, and these inconsistencies are the main problem with this\nresponse.\nAs already mentioned by Reviewer #2 (first comment on the first submission), in a project of\nthis magnitude it is inevitable that the material will be analysed by different groups of\nauthors and that the results will also be described and interpreted by different groups of\nauthors. Therefore, Reviewer #2 asks the corresponding authors to go through the entire\nmanuscript in order to make it more consistent.\nIn my opinion, this has not yet been done. I am even more disturbed by the fact that not\npossibly all analyses were carried out with the same material. Just a few examples:\n1) Supplementary Table 1 lists 514 samples, but at the beginning of the Results section, lines\n139-145, the authors mentioned 517 samples (297+198+22=517).\n2) Eight Banteng are mentioned in the main text, but only five are listed in Table S1.\n3) Table S11 lists eight banteng, but only B1, B2 and B3 are identical to Table S1. Banteng\nlike banteng08, banteng09, banteng_ypt2224, banteng_ypt2225, banteng_ypt2230 are not\ndescribed in the main material (Table S1).\n4) Heading of the important table S11 is wrong (source group twice).\n5) Cross-reference to Table S11 is wrong (see Supplementary Note 4 line 365).\n6) Supplementary Note 4 mentions gaur04, but there is no gaur04 in Table S1.\n7) Supplementary Note 4 mentions 40 core and pure Indicine cattle, but in Table S1 (last\ncolumn) we find only 20 Indicine cattle. Table S11 mentions 17 pure SAI animals. Most\nprobably we should add Sha3b, Thar1 and Har03 to this list and thus get 20 again and not\n40. The Response to my first comment also mentions 40 Core and pure Indicine cattle.\n8) Table S1 and Figure S7 do not use the same sorting. Therefore, the reader cannot assign\nthe results of the Admixture analysis to a particular Biosample ID.\nThese are just inconsistencies I spotted because I was particularly interested in the use of\npure reference panels in the introgression analyses. I am not in a position, and it is not my\njob, to check all the cross-references and all the figures, but from the above I get the\nimpression that the authors do not consider Reviewer #2's first comment. The authors do\nnot address it either.\nWith respect to ADMIXTURE, TreeMix, D- and f3-statistics.\nAll these analyses are performed at the genome or even population level, which does not\nnecessarily indicate the purity of the animal at the locus level. The analyses that check for\npossible admixture at the locus/segment level are more appropriate. This is particularly\nimportant for the Banteng and Gaur reference population.\nAs mentioned in this paper (https://www.nature.com/articles/s41467-018-05257-7 ),\nfactors like large difference in the number of individuals sampled per breed, genetic drift in\nthe population can affect its ancestry composition. Similarly, post-admixture drift can also\ndistort the signals of f3-statistics; this is also mentioned by the authors themselves in the\nlegend of Table S11. Further, in the description of f3-statistics (supplementary note 4), the\nauthors incorrectly mentioned the threshold of Z-score as 3; for the test to be significant the\nZ-score of f3-tests should be less than -3. Next, f3-test is better suited to group of animals\n(increases the power of test) rather than performing on each individual separately. I would\nalso recommend the authors to perform f3-tests on each SAI cattle breed as target with\nBanteng as P1 (one of the source population) and Zebu as P2 (another source population) to\ncheck the sensitivity of this test.\nMethods used to identify signature of selection such as composite likelihood ratio (CLR) and\niHS has increased power if the information about derived alleles is used. Since the authors\nhave already aligned Gaur, Banteng, Bison, Wisent, Yak and Buffalo to the cattle assembly,\nthey could have easily ascertain the ancestral alleles. Therefore, I am curious as to why the\nauthors decided not to use ancestral allele information for this analysis. Further, the authors\nhave used the results of Tajima’s D to support their signature of selection analysis, however,\nin some cases (for instance, Fig 2 and 3), it is hardly negative. Usually, the significance of\nTajima’s D, whether or not the obtained value significantly deviates from 0, is obtained\nusing beta distribution (table 2, of original Tajima’s D paper). Therefore, I suggest that the\nauthors either remove the results of Tajima’s D or justify it.\nFurther, parameters used in D-statistics are not described in sufficient details. Please include\nit in Table S3 too.\nThere are 84 introgressed and then positively selected genes in 42 candidate regions (line\n274, U50). These are validated by phylogenetic analysis. Supplementary Figures 20-26 show\nthe phylogenetic analysis of 27 of these regions. All these phylogenetic trees indicate\ncomplete fixation of the introgressed segment. The authors show that EAI cattle are\ndescended from pure SAI cattle, but in the 27 regions presented here, not a single EAI\nhaplotype clustered with its native group of origin. Such 100% turnover could be expected\nfor one or two strongly positively selected segments, but I would not expect this for all 27\nsegments shown in Figures S20-S26.\nIn what demographic scenario would such a 100% turnover be expected for 27 segments\npresented?\nHow strong must positive selection be to observe the phylogenetic trees shown in Figures\nS20-S26?\nIn all 27 trees, EUT and SAI haplotypes form monophyletic groups, and these are\nparaphyletic with respect to EAI, and there is not a single EAI haplotype that cluster with its\norigin group SAI. Could the Funder effect cause the observed patterns?\nIn the main text (lines 279-288), the positively selected segment in the proximal region of\nBTA25 is presented and discussed (fig. S26 and fig. 4). According to this, EAI animals in most\nof the EAI subpopulations studied carry the native SAI haplotype with a frequency of >0.5.\nThe haplotypes of presumed Banteng and Gaur origin show only an increased frequency. In\nmy opinion, this strongly contradicts the phylogeny shown in Supplementary Figure 26,\nwhere no EAI haplotype was clustered with SAI.\nIn addition to Banteng and Gaur, the authors implicate Kouprey, Gayal (data origin not\ngiven) and some phantom species as explanations for the observed diversity in the proximal\nsegment of BTA25. The authors should at least try to discuss an explanation other than the\nimplausible and complicated multi-species introgression and subsequent positive selection?\nIf I remember correctly, the globin gene family, and in particular the α-globin chain, is\nknown for a complicated evolution of paralogous genes, including recurrent mutations,\ngene conversion, co-evolution of different family members, and so on. Authors should check\nif complicated evolution of paralogues gene families could affect their conclusions.\nFigure S24 shows three chromosomal segments on BTA18. All three are located in the most\ncomplex region of the Bos taurusa and Bos indicus genomes (see\nhttps://doi.org/10.3168/jemandes2021-21625). The segmental duplication behaves\nsimilarly to the gene families mentioned above, and the region on BTA18 is very difficult to\nreconstruct correctly. The authors should check whether segmental duplication (in general)\ncould influence their conclusions.\nOther minor comments\nLine 76-77, the authors write …“was a direct consequence of a small number of genomic\nregions.” While it is conceivable that introgressed segments in EAI cattle did provide\nadaptive advantage, this statement is too strong especially considering that authors\nthemselves have identified genomic segments under selection that were not introgressed\nfrom its wild relatives. Selection based on standing genetic variation should at least be\ndiscussed as a possibility. Further, authors should consider removing the phrase “direct\nconsequence” as the authors did not provide any evidences showing direct and strong\ncausality between introgressed segments and adaptation in EAI cattle.\nLine 113,”..where have large local population are found” This sentence appears to be\ngrammatically incorrect, please rephrase it.\nLine 205, section, “Ancestral environmental adaptation of…”  The authors should at least\nprovide a sentence or two to give readers some information about the methods applied to\nidentify signature of selection in this section. I am emphasizing this because when I started\nreading this section I was under the impression that the same set of methods were applied\nhere as applied in the previous section, i.e. Fst, theta and XP-EHH, but it was not until the\nend of this section that it was mention that a different set of methods were applied, i.e. CLR\nand iHS. Moreover, the naming of the sections are also confusing, “the ancestral adaptation\nof indicine cattle”, “Ancestral environmental adaptation of south Asian indicine cattle”, “the\nindicine adaptation to the tropical, humid environment”. Ideally, signature of selection\nanalysis should be performed without a priori expectation because it is very difficult to\nconclusively link the identified genes under selection with the breed characteristics; usually,\nit is done by functional genomics and by studying the literature extensively. Therefore, such\nnaming of the sections is suitable to discussion but not to result section.\nReviewer #2 (Remarks to the Author):\nOverall, the manuscript has improved significantly. I appreciate the author’s willingness to\nmake changes and incorporate suggestions. There are three documents as I see it, the main\nmanuscript, the supplement and the response to previous review, all of which are lengthy. I\nhave comments for each, which I’ll address individually, but the most pressing issue relates\nto the ROH analysis under Supplementary Note 2: Genetic diversity. I’ll start there because I\nthink this needs addressed.\nResponse L358-400 and supplement L273-280. I still do not believe your ROH analysis is\naccurate. Even with your new analysis, you show a large number of individuals with a total\ngenome length contained within ROH to be >1 Gb, or over a third of the genome. That is\nsimply not realistic. There has been a large amount of prior literature on the subject, largely\nfrom SNP-chip data, and all of them show the “extreme” samples having under 1 Gb in ROH.\nI believe this [1] is an informative publication.\n[1] https://bmcgenomdata.biomedcentral.com/articles/10.1186/1471-2156-13-70\nI believe Figure 5 from [1] is particularly useful in that it shows the same trend you are\nseeing in your data but the scale is quite different. This is particularly important because if\nthe current manuscript is published it will present conflicting values because I believe you\nare including far too much of the genome in ROH. It is generally accepted that these types of\nROH analysis are measuring recent inbreeding, which you acknowledge on L160 of the\nmanuscript, where the size of ROH regions is proportional to the inbreeding coefficient. As\nthe size of a ROH decreases you are estimating IBD further back in time. This may be\nappropriate in some cases, even perhaps the current manuscript, but your ability to\ndifferentiate old IBD regions from random chance or other population forces requires\nincreased samples sizes, which you do not have. As I said, on it’s face, I do not believe your\nROH results. The question is, what is causing them to be elevated? In order to try to figure\nthis out, I extracted 187 samples from your study that are also present within my UMAG1\ncall set, which is based on ~5500 genomes and VQSR. For these 187 samples, I extracted two\nsets of variants, A) all bi-allelic SNPS that passed VQSR and B) the positions contained on the\nIllumina Bovine HD and GGPF250 genotyping assays in order to compare with the results\nfrom [1].\nA) “bcftools view --threads $CPU -f PASS -m2 -M2 -v snps -S $InputSamples --force-samples -\nOb -o $c.$Prefix.bcf $InputBcf”\nB) “bcftools view --threads $CPU -R $c.variants2get.list --trim-alt-alleles -Ob -o $c.$Prefix.bcf\n$InputBcf”\nDataset A produced 125M variants for these 187 samples.\nbcftools +counts ChenAll.bcf\nNumber of samples: 187\nNumber of SNPs: 125813861\nNumber of INDELs: 0\nNumber of MNPs: 0\nNumber of others: 0\nNumber of sites: 125813861\nDataset B produced ~850K variants.\nbcftools view -S chen_sample.list --force-samples ChenChip.bcf -Ou | bcftools +counts\nNumber of samples: 187\nNumber of SNPs: 850772\nNumber of INDELs: 64337\nNumber of MNPs: 0\nNumber of others: 0\nNumber of sites: 897036\nThe bcf files were then converted to plink bed/bim/bam:\nplink --bcf ChenAll.bcf --const-fid 0 --chr-set 29 --allow-extra-chr --out ChenAll\nplink --bcf ChenChip.bcf --const-fid 0 --chr-set 29 --allow-extra-chr --out ChenChip\nFinally, the ROH analysis was performed using various parameter values.\nThe first run used the values specified on L279-280 of main paper and L389-396 of\nResponse. I will refer to this as “ChenParameters”.\nplink --bfile ChenAll --chr-set 29 --chr 1-29 --homozyg --homozyg-gap 1000 --homozyg-kb\n100 --homozyg-snp 50 --homozyg-window-het 3 --homozyg-window-snp 50 --homozyg-\nwindow-threshold 0.05 --out ChenAll\nThe second run was the same as above but changed --homozyg-window-het 3 to 1. I refer to\nthis as Chen1het.\nThe third run was the same as above but changed --homozyg-window-het 3 --> 1, --\nhomozyg-snp 50 --> 200 --homozyg-window-snp 50 --> 100 --homozyg-snp 50 --> 200. I refer\nto this as ChenStrict.\nThe loci from the SNP-chips were run the same as the second analysis above which used the\nauthor’s parameters but changed --homozyg-window-het 3 --> 1. I refer to this as Chip1het.\nAll of the results from the *.hom.indiv plink files were copied into Excel and visualized. See\nfile ChenROHsubmitted.xlsx which should be available with this review. It is clear from\nFigure 1 that I was able to recreate the trend that is shown in Supp. Figures 2 & 3 with my\n“ChenParameters” in blue. As I contend, these parameters still include too many ROH and\nencompass too much of the genome per individual. The “Chen1het” in orange only change -\n-homozyg-window-het from 3 to 1 which produces values that more closely resemble prior\nliterature and what we know about these populations. Therefore, I believe the root problem\nwith the analyses, as presented, is this parameter. In hindsight, this makes sense as this\nparameter is meant to allow for genotyping error. By setting this parameter to 3 and\nrequiring only 100 loci within a window, I believe you are effectively saying that you have\nlittle confidence in your genotype calls because you are allowing 3 errors (3% error rate) to\nstill call a region as homozygous. If you are confident in your genotype calls, and new\nmutations are negligible (which I think is true here) then a value of 1 is more appropriate\nand in fact yields results closer to expectations. In order to more fully explore this, I also\ncreated a more strict dataset (ChenStrict) which required more loci while still maintaining 1\nhet position (grey). In order to provide a direct comparison to snp-chip data I also present\nthe same as “Chip1het” in red. As you can see, the “Chip1het” is comparable to prior\nliterature from using these assays. My final interpretation of these data is that the\nparameters used in the revision are still not appropriate, the main parameter to change is --\nhomozyg-window-het and the correct values for the other parameters are somewhere\nbetween what was used and the “strict” values that I used. This is a fairly trivial analysis that\nthe authors should perform and update the manuscript.\nAs before, I will provide comments on the Supplement by referring to SL as supplement line\nnumber.\nSL261-264: “The whole genome sequencing data from six other bovine species were\nmapped in the same way. We used the 67,162,108 SNPs as a reference list to genotype the\ncombined set of 495 cattle samples and 22 samples of six other bovine species, resulting in\n67,145,163 SNP data with wild species data.”\nIt is unclear to me what you mean by “…as a referencelist…”? My initial interpretation of\nthis is that you added additional samples, called genotypes on those additional samples,\nand then just extracted the 67M positions from the prior variant call set. Is this correct? The\nproper way to perform this analysis would be to recall genotypes from *all* of your samples\nstarting with the combined GVCFs from SL253. This confusion is exacerbated by the third\ndescription in Supp. Table 3 which says “A total of 64,475,272 SNPs called from all 495 cattle\nand 22 genomes from other six bovine species.” Please clarify exactly what you did here.\nSL266 Supplementary Table 3 This table is very informative and helpful to the reader and\nthis reviewer. However, now that I have this table, it presents new problems. For the second\nline of this table where it says “Imputed data of…” I assume this represents the phasing and\nimputation that you did with Beagle, correct? If so, please add the word “phased” to this\ndescription. This is a very important point later.\nNone of these numbers seem to add up, going back to my comment directly above. It seems\nto me that the proper way to analyze all of these data is to do joint genotyping on all samples\nsimultaneously, which would be your 514 samples presented in Supp. Table 1.\n(Note that the second line of this table, and in the manuscript and supplement you refer to\n495+22 which is 517 but you only have 514 in Table1.) This will result in a total number of\nvariants N. All analyses performed after this will involve some level of filtering or selection\nof loci but they all will go back to this original N variants. As it stands, this table helped clear\nup some issues but it presented more. For example the 6th description in the number of\nsamples column, “All 495 cattle and one yak”. Table 1 lists 3 yak while SL237 says two yak.\nTable 1 has 5 banteng while SL237 says eight. Figure 4 (F) of the main manuscript mentions\nKouprey but I do not see any in table 1 as Kouprey or Bos sauveli. Where did this sample\ncome from? Again, many of these numbers are not consistent.\nSL285 and 300 The PCA and admixture analysis used -indep-pairwise 50 10 0.1 while the NJ\nand ML phylogeny used --indep-pairwise 50 5 0.1. Why would you use different parameters\nfor two different analyses to achieve the same objective of minimizing LD between the\nvariants you use? By using different parameter sets you effectively make a direct\ncomparison of these two analyses impossible because you changed the loci that you are\nusing.\nSL312 “…XP-EHH score using selscan v.1.1 with default settings.” Selscan software did not\nhave the ability to use unphased data until v2.0. Therefore, I assume that you are relying on\nthe Beagle phased and imputed data for these analyses, is that correct? If so, this is why it is\nimportant to add the word phased to Supp. Table 1 because without doing the phasing in\nBeagle you could not have used selscan v.1.1.\nSL329 “… using SweepFinder2.”\nFrom the manual, it appears that this requires a recombination rate file. What did you use\nas the recombination rate between loci? My assumption is that you assumed a constant\nrecombination rate of 1 cM/Mb since you do not actually have a recombination map. If that\nis the case, please state that in the manuscript. This software also requires a “B-value” for\neach variant. Where did you get these values from?\nSL354: Why were EAI grouped into a single population without looking into recent specific\nbreed introgression as the potential source of admixed genotypes? Often times the\ndirection of admixture events is difficult to estimate. Is the signal of introgression the result\nof the modern banteng/gaur samples having some amount of domestic cattle introgression\nin their genomes?\nSL360-S365: The f3 data does not appear to be presented in Supplementary table 10, There\nappears to be a numbering issue of tables in the supplementary data document, with a shift\nof +1 after supplementary table 7.\nMain manuscript referring to ML as manuscript line number:\nML 76 “…was a direct consequence of…” That’s a fairly strong claim. I think you evidence\nconsistent with a hypothesis but I question whether it represents a direct consequence.\nML113 remove the word have.\nML261 Supp. Figure 16 has multiple pages of figures with the same panel labels but only 1\nlegend. Please have a look at this as it appears it may have been an error when preparing\nthe files for submission.\nML264-266 “To identify regions in the EAI genomes that were likely under selection, we\nused a statistic to detect positively selected and introgressed genes (PSIGs)” Are these\nactually detecting genes or regions? You use both terms in the same sentence.\nML351 Supp. Fig. 38 This is a much better figure than the previous version!\nML512 “…gene transfer format (GTF) (GCF_002263795.1) file.” This is an accession for the\ngenome assembly, it is NOT an accession for a gene set. I assume that you either used\nannotations from NCBI or Ensembl, both of which have many different versions which are\nregularly updated. You need to specify which annotations you used (NCBI or Ensembl) and\nexactly which version of the annotation.\n2 [Comment of Reviewer #1]\n5 from pure Indicine cattle originating in the South Asian core area (i.e. SAI, Indus Valley).\n7 successful adaptation to hot and humid environments occurred through introgression of\n9 environment is estimated to be about 3000 years, i.e. 600 generations of cattle. That mean the\n10 first contact to and the introgression of Banteng and/or Gaur should be about 600 generations\n11 ago or younger. These results are supported by cumulative evidence from different analyses.\n13 becoming more popular and many signals are interpreted as introgression, often without an\n14 alternative hypothesis.\n16 supported by the GO and KEGG enrichment results but would of course need additional\n17 evidence from functional studies.\n18 [Comment of Reviewer #1]\n20 important for the above mentioned main message of this manuscript. I asked questions like:\n21 1) Do the authors ensure that the SAIs are pure indicines without a genomic taurine segment?\n22 2) Do the authors ensure that the Banteng samples are pure and do not contain a genomic\n23 segment from indicine or taurine cattle?\n24 3) Do the authors ensure that the Gaur samples are pure and do not contain a genomic segment\n25 from indicine or taurine cattle?\n27 inconsistencies that worry me, and these inconsistencies are the main problem with this\n28 response.\n30 a group of proximate cross-fertile species. However, according to your suggestions, we have\n31 checked the D statistics among indicine cattle, taurine cattle, banteng and gaur samples to\n32 select pure SAI samples without taurine, banteng or gaur introgression, select pure taurine\n33 cattle without SAI cattle, banteng or gaur introgression, and select banteng or gaur without SAI\n34 and taurine introgression. Finally, we selected 15 pure SAI cattle, 15 taurine cattle, 4 banteng,\n35 and 2 gaur and used them for introgression analysis (Supplementary Table 10 and\n38 introgression analysis. For SAI cattle, we used the three tree topologies of D (SAI individual,\n40 buffalo), and D (SAI individual, SAI individual, gaur, buffalo) to select the SAI samples without\n41 any gene flow from taurine cattle, banteng or gaur. For taurine cattle, we used three tree\n42 topologies of D (taurine individual, taurine individual; SAI, buffalo), D (taurine individual,\n43 taurine individual; banteng individual, buffalo), and D (taurine individual, taurine individual;\n44 gaur individual, buffalo) to select taurine samples without any gene flow from SAI cattle,\n45 banteng or gaur. For banteng, we used two tree topologies of D (banteng individual, banteng\n46 individual; SAI individual, buffalo) and D (banteng individual, banteng individual; taurine\n47 individual, buffalo) to select banteng samples without any gene flow from taurine or SAI cattle.\n49 buffalo) and D (gaur individual, gaur individual; taurine individual, buffalo) to select gaur\n50 samples without any gene flow from taurine or SAI cattle. We used all combinations of\n51 individuals to calculate the D statistic. We finally selected a panel of 15 pure SAI cattle, 15\n52 taurine cattle, 4 banteng, and 2 gaur samples with a |Z score| < 3 for RFmix analysis, D statistic,\n55 [Comment of Reviewer #1]\n57 this magnitude it is inevitable that the material will be analysed by different groups of authors\n58 and that the results will also be described and interpreted by different groups of authors.\n60 order to make it more consistent.\n62 all analyses were carried out with the same material. Just a few examples:\n63 1) Supplementary Table 1 lists 514 samples, but at the beginning of the Results section, lines\n64 139-145, the authors mentioned 517 samples (297+198+22=517).\n66 495 cattle genomes and 22 wild bovine species genomes. We also used one gayal and one\n67 ancient kouprey sample for the introgression region analysis of BTA25. We added this\n68 information in the Supplementary Note 4.\n69 “For the analysis of the introgressed region of BTA25 (0.21-0.26 Mb), we also used a gayal\n70 sample and an ancient kouprey sample to detect its origin. The coverage of gayal and kouprey\n71 was 17.32× and 1.4×, respectively. Due to the hybrid origin of gayal and low coverage of the\n72 kouprey genome, we did not examine possible introgression of gayal and kouprey. The publicly\n73 available sequences were downloaded from China National GeneBank (CNGB) and the SRA\n74 with the following project accession numbers: CRX165997 (gayal, YD4) and PRJNA764746\n75 (kouprey).” Considering the low coverage of the kouprey genome and the hybrid ancestry of\n76 gayal, we did not further analyze the gene flow between these two species and EAI cattle.\n77 [Comment of Reviewer #1]\n78 2) Eight Banteng are mentioned in the main text, but only five are listed in Table S1.\n80 dataset. According to your comment, we updated Supplementary Table S1 and added all\n81 information on the 8 banteng, and we only used 4 banteng and gaur samples in the introgression\n82 analysis.\n83 [Comment of Reviewer #1]\n84 3) Table S11 lists eight banteng, but only B1, B2 and B3 are identical to Table S1. Banteng like\n85 banteng08, banteng09, banteng_ypt2224, banteng_ypt2225, banteng_ypt2230 are not\n86 described in the main material (Table S1).\n88 our samples to ensure that the names are consistent in all materials.\n89 [Comment of Reviewer #1]\n90 4) Heading of the important table S11 is wrong (source group twice).\n92 and we selected pure SAI cattle, taurine cattle, banteng, and gaur for f analysis.\n3\n93 [Comment of Reviewer #1]\n94 5) Cross-reference to Table S11 is wrong (see Supplementary Note 4 line 365).\n96 [Comment of Reviewer #1]\n97 6) Supplementary Note 4 mentions gaur04, but there is no gaur04 in Table S1.\n99 [Comment of Reviewer #1]\n100 7) Supplementary Note 4 mentions 40 core and pure Indicine cattle, but in Table S1 (last column)\n101 we find only 20 Indicine cattle. Table S11 mentions 17 pure SAI animals. Most\n102 probably we should add Sha3b, Thar1 and Har03 to this list and thus get 20 again and not\n103 40. The Response to my first comment also mentions 40 Core and pure Indicine cattle.\n105 cattle, 4 banteng, and 2 gaur for introgression analysis. According to your suggestions, we have\n106 checked the D statistic among indicine, taurine, banteng and gaur samples. We have corrected\n107 our mistakes in Supplementary Note 4 and Supplementary Tables 10 and 11.\n108 [Comment of Reviewer #1]\n109 8) Table S1 and Figure S7 do not use the same sorting. Therefore, the reader cannot assign the\n110 results of the Admixture analysis to a particular Biosample ID.\n112 and Supplementary Figure 7 to ensure consistency.\n114 pure reference panels in the introgression analyses. I am not in a position, and it is not my\n115 job, to check all the cross-references and all the figures, but from the above I get the\n116 impression that the authors do not consider Reviewer #2's first comment. The authors do\n117 not address it either.\n120 [Comment of Reviewer #1]\n122 at the genome or even population level, which does not necessarily indicate the purity of the\n123 animal at the locus level. The analyses that check for possible admixture at the locus/segment\n124 level are more appropriate. This is particularly important for the Banteng and Gaur reference\n125 population.\n127 groups, we calculated the D statistic and selected 15 pure South Asian indicine (SAI) cattle, 15\n128 taurine cattle, 4 banteng, and 2 gaur for introgression analysis. We agree with your point, and\n129 our results in Supplementary Figs. 20-25 show locus-specific trees, which show if a fragment\n130 introgressed from, for example, banteng has a non-banteng origin.\n131 [Comment of Reviewer #1]\n133 large difference in the number of individuals sampled per breed, genetic drift in the population\n134 can affect its ancestry composition. Similarly, post-admixture drift can also distort the signals\n135 of f3-statistics; this is also mentioned by the authors themselves in the legend of Table S11.\n137 mentioned the threshold of Z-score as 3; for the test to be significant the Z-score of f3-tests\n138 should be less than -3. Next, f3-test is better suited to group of animals (increases the power of\n139 test) rather than performing on each individual separately. I would also recommend the authors\n140 to perform f3-tests on each SAI cattle breed as target with Banteng as P1 (one of the source\n141 population) and Zebu as P2 (another source population) to check the sensitivity of this test.\n143 interpreted with caution. However, our Admixture pattern reproduce the divergence of\n144 geographically separated clusters: European taurine, African taurine, African indicine, South\n146 groups, we calculated the D statistic and selected 15 pure South Asian indicine (SAI) cattle, 15\n147 taurine cattle, 4 bangteng, and 2 gaur samples. We also calculated f statistics for each SAI\n3\n148 breed as targets with banteng, gaur, taurine, and other SAI breeds. For the f statistic, if the Z\n3\n149 score (Z ≤−3.0) is significantly negative, test population C has admixture from both reference\n3\n150 populations A and B. All f statistics were positive, indicating that there was no evidence of\n3\n151 admixture (Supplementary Table 11).\n152 [Comment of Reviewer #1]\n154 iHS has increased power if the information about derived alleles is used. Since the authors have\n155 already aligned Gaur, Banteng, Bison, Wisent, Yak and Buffalo to the cattle assembly, they\n156 could have easily ascertain the ancestral alleles. Therefore, I am curious as to why the authors\n157 decided not to use ancestral allele information for this analysis.\n160 and derived allele state is needed for each SNP. In our analysis, the ancestral allele was defined\n161 as the allele fixed in the swamp buffalo that was included in the genotype call set, and the\n162 ambiguous SNP was discarded. We also updated the iHS and CLR results in our manuscript.\n164 [Comment of Reviewer #1]\n166 analysis, however, in some cases (for instance, Fig 2 and 3), it is hardly negative. Usually, the\n167 significance of Tajima’s D, whether or not the obtained value significantly deviates from 0, is\n168 obtained using beta distribution (table 2, of original Tajima’s D paper). Therefore, I suggest that\n169 the authors either remove the results of Tajima’s D or justify it.\n171 in Table S3 too.\n173 parameters used in calculating the D statistic in Supplementary Table 3.\n174 [Comment of Reviewer #1]\n176 274, U50). These are validated by phylogenetic analysis. Supplementary Figures 20-26 show\n177 the phylogenetic analysis of 27 of these regions. All these phylogenetic trees indicate\n178 complete fixation of the introgressed segment. The authors show that EAI cattle are descended\n179 from pure SAI cattle, but in the 27 regions presented here, not a single EAI haplotype clustered\n180 with its native group of origin. Such 100% turnover could be expected for one or two strongly\n181 positively selected segments, but I would not expect this for all 27 segments shown in Figures\n184 presented?\n188 paraphyletic with respect to EAI, and there is not a single EAI haplotype that cluster with\n189 its origin group SAI. Could the Funder effect cause the observed patterns?\n191 our previous revision, we used the same color for indicine or taurine haplotypes. In the indicine\n192 group, we did find EAI haplotypes clustered with SAI, but we did not color them.\n194 introgressed regions, and used 5 SAI cattle, 5 taurine cattle, 4 banteng, and 2 gaur samples and\n195 other wild species samples to construct phylogenetic trees and validate the introgression.\n197 introgressed haplotypes of EAI cattle were clustered with banteng or gaur, while no introgressed\n198 haplotypes were clustered with South Asian indicine cattle (Supplementary Figs. 20-25).\n199 [Comment of Reviewer #1]\n202 of the EAI subpopulations studied carry the native SAI haplotype with a frequency of >0.5.\n204 opinion, this strongly contradicts the phylogeny shown in Supplementary Figure 26, where no\n208 haplotypes, and new figures showed that non-introgressed haplotypes of EAI cattle were\n209 clustered with SAI haplotypes. In Supplementary Figs. 20-25, we first specified the length of\n210 the introgressed regions and used 5 pure SAI cattle, 5 taurine cattle, 4 banteng, 2 gaur, 3 yak,\n211 2 bison, 2 wisent, 2 buffalo, and 80 EAI samples to construct phylogenetic trees. Our results\n212 showed that introgressed haplotypes of EAI were clustered with banteng or gaur, while no\n213 introgrssed haplotypes were clustered with SAI haplotypes. We also found that some EAI\n214 haplotypes also clustered with taurine cattle, which showed that some EAI also have taurine\n215 ancestry.\n217 wisent-yak outgroup or gaur-banteng is variable. In some trees, banteng and gaur are\n218 separated. We also checked that the genomic segments were demarcated accurately and\n219 contained only the introgressed fragments. Please see Supplementary Figs. 20-25.\n220 [Comment of Reviewer #1]\n222 and some phantom species as explanations for the observed diversity in the proximal segment\n223 of BTA25. The authors should at least try to discuss an explanation other than the\n224 implausible and complicated multi-species introgression and subsequent positive selection?\n226 everywhere the species share territory. The cluster patterns indeed suggest sequences not\n227 originating from zebu, banteng, gaur or kouprey, but this may reflect, especially for the gaur,\n228 their limited sampling. For the selected genes, we have revised our manuscript and now refer\n229 to the selected genes as frequently introgressed genes.\n230 [Comment of Reviewer #1]\n232 known for a complicated evolution of paralogous genes, including recurrent mutations,\n233 gene conversion, co-evolution of different family members, and so on. Authors should check\n234 if complicated evolution of paralogues gene families could affect their conclusions.\n236 shows two HBA genes that encode identical HbA subunits. We did not find reports on HbA\n237 rearrangements or paralogs that would invalidate our results.\n238 [Comment of Reviewer #1]\n240 complex region of the Bos taurusa and Bos indicus genomes (see\n241 https://doi.org/10.3168/jemandes2021-21625). The segmental duplication behaves similarly to\n242 the gene families mentioned above, and the region on BTA18 is very difficult to reconstruct\n243 correctly. The authors should check whether segmental duplication (in general)\n244 could influence their conclusions.\n246 https://doi.org/10.3168/jds.2021-21625. This paper indeed mentions SVs near a QTL region on\n247 (BTA18, 57,816,000-59,430,000), close to the regions shown in Supplementary Fig. 24 (now in\n249 we are studying, but one of their samples has a deletion of 59123315 to 61313922. This would\n250 delete all segments on which the three BTA18 trees in Supplementary Fig. 23 are based, but\n251 does not seem to influence the phylogenetic identification of introgression in other samples. In\n252 the legend of Fig. S23 we have mentioned the occurrence of this SV with reference to Dachs et\n253 al (2023).\n254 [Comment of Reviewer #1]\n257 regions.” While it is conceivable that introgressed segments in EAI cattle did provide\n258 adaptive advantage, this statement is too strong especially considering that authors\n259 themselves have identified genomic segments under selection that were not introgressed\n260 from its wild relatives. Selection based on standing genetic variation should at least be\n261 discussed as a possibility. Further, authors should consider removing the phrase “direct\n262 consequence” as the authors did not provide any evidences showing direct and strong\n263 causality between introgressed segments and adaptation in EAI cattle.\n265 [Comment of Reviewer #1]\n267 grammatically incorrect, please rephrase it.\n269 [Comment of Reviewer #1]\n271 provide a sentence or two to give readers some information about the methods applied to\n272 identify signature of selection in this section. I am emphasizing this because when I started\n273 reading this section I was under the impression that the same set of methods were applied here\n274 as applied in the previous section, i.e. Fst, theta and XP-EHH, but it was not until the end of\n275 this section that it was mention that a different set of methods were applied, i.e. CLR\n276 and iHS. Response: Thank you for this suggestion. We have added the description in the first and last\n277 paragraphs.\n279 “Throughout the history of migration and admixture of indicine cattle, genomic regions under\n280 selection might have been lost in specific indicine groups. We therefore also performed a test\n281 for positive selection signatures in SAI cattle using θπ, iHS, CLR, and F estimates based on\nST\n282 the comparison between SAI and non-SAI groups.”\n284 “Additionally, we compared the light- and dark-coated SAI breeds, such as the white-coated\nST\n286 shared selective sweeps around pigmentation loci, e.g., LEF1 and ASIP, in the light-coated\n287 indicine breeds (Fig. 2). This selection pressure may have been favored or driven by high\n288 temperatures and intense solar radiation and/or human preferences. Across the whole genomes,\n289 the CLR and iHS analyses revealed 368 regions overlapping with 477 genes present in AFI and\n291 equally important for AFI cattle”\n293 cattle”, “Ancestral environmental adaptation of south Asian indicine cattle”, “the\n294 indicine adaptation to the tropical, humid environment”. Ideally, signature of selection\n295 analysis should be performed without a priori expectation because it is very difficult to\n296 conclusively link the identified genes under selection with the breed characteristics; usually,\n297 it is done by functional genomics and by studying the literature extensively. Therefore, such\n298 naming of the sections is suitable to discussion but not to result section.\n300\n302 [Comment of Reviewer #2]\n304 make changes and incorporate suggestions. There are three documents as I see it, the main\n305 manuscript, the supplement and the response to previous review, all of which are lengthy. I have\n306 comments for each, which I’ll address individually, but the most pressing issue relates to the\n308 this needs addressed.\n310 accurate. Even with your new analysis, you show a large number of individuals with a total\n311 genome length contained within ROH to be >1 Gb, or over a third of the genome. That is simply\n312 not realistic. There has been a large amount of prior literature on the subject, largely from SNP-\n313 chip data, and all of them show the “extreme” samples having under 1 Gb in ROH. I believe\n314 this [1] is an informative publication.\n315 [1]https://bmcgenomdata.biomedcentral.com/articles/10.1186/1471-2156-13-70 I believe\n317 data but the scale is quite different. This is particularly important because if the current\n318 manuscript is published it will present conflicting values because I believe you are including\n319 far too much of the genome in ROH. It is generally accepted that these types of ROH analysis\n320 are measuring recent inbreeding, which you acknowledge on L160 of the manuscript, where\n321 the size of ROH regions is proportional to the inbreeding coefficient. As the size of a ROH\n322 decreases you are estimating IBD further back in time. This may be appropriate in some cases,\n323 even perhaps the current manuscript, but your ability to differentiate old IBD regions from\n324 random chance or other population forces requires increased samples sizes, which you do not\n325 have. As I said, on it’s face, I do not believe your ROH results. The question is, what is causing\n326 them to be elevated? In order to try to figure this out, I extracted 187 samples from your study\n327 that are also present within my UMAG1 call set, which is based on ~5500 genomes and VQSR.\n329 and B) the positions contained on the Illumina Bovine HD and GGPF250 genotyping assays in\n330 order to compare with the results from [1].\n332 -Ob -o $c.$Prefix.bcf $InputBcf”\n334 $InputBcf”\n336 bcftools +counts ChenAll.bcf\n344 bcftools view -S chen_sample.list --force-samples ChenChip.bcf -Ou | bcftools +counts\n352 plink --bcf ChenAll.bcf --const-fid 0 --chr-set 29 --allow-extra-chr --out ChenAll\n353 plink --bcf ChenChip.bcf --const-fid 0 --chr-set 29 --allow-extra-chr --out ChenChip\n357 plink --bfile ChenAll --chr-set 29 --chr 1-29 --homozyg --homozyg-gap 1000 --homozyg-kb\n358 100 --homozyg-snp 50 --homozyg-window-het 3 --homozyg-window-snp 50 --homozyg-\n359 window-threshold 0.05 --out ChenAll\n361 this as Chen1het.\n363 snp 50 --> 200 --homozyg-window-snp 50 --> 100 --homozyg-snp 50 --> 200. I refer to this as\n366 author’s parameters but changed --homozyg-window-het 3 --> 1. I refer to this as Chip1het.\n368 file ChenROHsubmitted.xlsx which should be available with this review. It is clear from Figure\n369 1 that I was able to recreate the trend that is shown in Supp. Figures 2 & 3 with my\n370 “ChenParameters” in blue. As I contend, these parameters still include too many ROH and\n371 encompass too much of the genome per individual. The “Chen1het” in orange only change --\n372 homozyg-window-het from 3 to 1 which produces values that more closely resemble prior\n373 literature and what we know about these populations. Therefore, I believe the root problem with\n374 the analyses, as presented, is this parameter. In hindsight, this makes sense as this parameter is\n375 meant to allow for genotyping error. By setting this parameter to 3 and requiring only 100 loci\n376 within a window, I believe you are effectively saying that you have little confidence in your\n377 genotype calls because you are allowing 3 errors (3% error rate) to still call a region as\n378 homozygous. If you are confident in your genotype calls, and new mutations are negligible\n379 (which I think is true here) then a value of 1 is more appropriate and in fact yields results closer\n380 to expectations. In order to more fully explore this, I also created a more strict dataset\n381 (ChenStrict) which required more loci while still maintaining 1 het position (grey). In order to\n382 provide a direct comparison to snp-chip data I also present the same as “Chip1het” in red. As\n383 you can see, the “Chip1het” is comparable to prior literature from using these assays. My final\n384 interpretation of these data is that the parameters used in the revision are still not appropriate,\n385 the main parameter to change is --homozyg-window-het and the correct values for the other\n386 parameters are somewhere between what was used and the “strict” values that I used. This is a\n387 fairly trivial analysis that the authors should perform and update the manuscript.\n389 your efforts! We agree that the “strict” values (ChenStrict) of the ROH detection method are\n390 more accurate, so we used the “strict” value parameters outlined in your comments to detect\n391 the ROH in our cattle data. The final parameters were set to a minimum length of 100 kb, a\n392 scanning window size of 100 SNPs, a minimum density threshold of 200 SNPs, a large gap of\n393 1000 kb, a maximum number of heterozygous SNPs in the scanning window of 1, and a scanning\n394 window threshold level of 0.05. These settings yielded expected number (maximum number was\n395 3259) and total length (maximum length was 1,138,710 Mb) of ROH (Supplementary Fig. 2).\n397\n399 geographic group.\n400 [Comment of Reviewer #2]\n402 number.\n404 in the same way. We used the 67,162,108 SNPs as a reference list to genotype the combined set\n405 of 495 cattle samples and 22 samples of six other bovine species, resulting in 67,145,163 SNP\n406 data with wild species data.”\n408 is that you added additional samples, called genotypes on those additional samples, and then\n409 just extracted the 67M positions from the prior variant call set. Is this correct? The proper way\n410 to perform this analysis would be to recall genotypes from *all* of your samples starting with\n411 the combined GVCFs from SL253. This confusion is exacerbated by the third description in\n413 from other six bovine species.” Please clarify exactly what you did here.\n415 sentence: “We genotyped the combined set of 495 cattle samples and 22 samples of six other\n416 bovine species, and then extracted the 67,162,108 SNPs. After filtering out the non-biallelic\n418 [Comment of Reviewer #2:]\n420 reviewer. However, now that I have this table, it presents new problems. For the second line of\n421 this table where it says “Imputed data of…” I assume this represents the phasing and\n422 imputation that you did with Beagle, correct? If so, please add the word “phased” to this\n423 description. This is a very important point later.\n425 imputed and phased data. We revised the text as follows: “We used phased and imputed SNP\n426 data”.\n428 me that the proper way to analyze all of these data is to do joint genotyping on all samples\n429 simultaneously, which would be your 514 samples presented in Supp. Table 1. (Note that the\n430 second line of this table, and in the manuscript and supplement you refer to 495+22 which is\n431 517 but you only have 514 in Table1.) This will result in a total number of variants N. All\n432 analyses performed after this will involve some level of filtering or selection of loci but they all\n433 will go back to this original N variants. As it stands, this table helped clear up some issues but\n434 it presented more. For example the 6th description in the number of samples column, “All 495\n435 cattle and one yak”. Table 1 lists 3 yak while SL237 says two yak. Table 1 has 5 banteng while\n437 in table 1 as Kouprey or Bos sauveli. Where did this sample come from? Again, many of these\n438 numbers are not consistent.\n441 revised the number of samples. For the 6th description, we selected only one yak for the outgroup,\n442 but to ensure consistency, we reanalyzed the data using three yak for the outgroup. We have\n443 eight banteng and five gaur in our dataset, and we selected four pure banteng and 2 gaur\n444 samples for the introgression analysi. We apologize for our mistakes in Supplementary Table 1\n445 and we updated it, and we also added the ancient kouprey and gayal information to\n447 “For the analysis of the introgressed region of BTA25 (0.21-0.26 Mb), we also used a gayal\n448 sample and an ancient kouprey sample to detect its origin. The coverage of gayal and kouprey\n449 was 17.32× and 1.4×, respectively. Due to the hybrid origin of gayal and low coverage of the\n450 kouprey genome, we did not examine possible introgression of gayal and kouprey. The publicly\n451 available sequences were downloaded from China National GeneBank (CNGB) and the SRA\n452 with the following project accession numbers: CRX165997 (gayal, YD4) and PRJNA764746\n453 (kouprey).”\n454 [Comment of Reviewer #2:]\n456 and ML phylogeny used --indep-pairwise 50 5 0.1. Why would you use different parameters\n457 for two different analyses to achieve the same objective of minimizing LD between the\n458 variants you use? By using different parameter sets you effectively make a direct\n459 comparison of these two analyses impossible because you changed the loci that you are\n460 using.\n462 10 0.1 for PCA and ML phylogenetic tree reconstruction, and all 67,162,108 autosomal SNPs\n463 were used to construct the NJ tree.\n464 [Comment of Reviewer #2]\n466 have the ability to use unphased data until v2.0. Therefore, I assume that you are relying on the\n468 important to add the word phased to Supp. Table 1 because without doing the phasing in Beagle\n469 you could not have used selscan v.1.1.\n471 software, and we have added this information to Supplementary Table 1.\n472 [Comment of Reviewer #2]\n475 recombination rate between loci? My assumption is that you assumed a constant recombination\n476 rate of 1 cM/Mb since you do not actually have a recombination map. If that is the case, please\n477 state that in the manuscript. This software also requires a “B-value” for each variant. Where did\n478 you get these values from?\n482 and the “recombination rate” and “B-value” are not required for each method. We applied\n483 “Scan for selective sweeps with pre-computed empirical spectrum” in SweepFinder2, and the\n484 parameter was “./SweepFinder2 --lu GridFile FreqFile SpectFile OutFile”.\n485 [Comment of Reviewer #2]\n487 breed introgression as the potential source of admixed genotypes?\n489 sum and added the D values. The EAI breeds have similar values.\n490 [Comment of Reviewer #2]\n492 introgression the result of the modern banteng/gaur samples having some amount of domestic\n493 cattle introgression in their genomes?\n495 detected by RFMix or by phylogenetic analysis.\n496 [Comment of Reviewer #2]\n498 appears to be a numbering issue of tables in the supplementary data document, with a shift of\n499 +1 after supplementary table 7.\n501 f3 results to Supplementary Tables 10 and 11, and the Supplementary Dataset1. We have double-\n502 checked the numbers of Supplementary tables and figures.\n503 [Comment of Reviewer #2]\n506 consistent with a hypothesis but I question whether it represents a direct consequence.\n508 [Comment of Reviewer #2]\n511 [Comment of Reviewer #2]\n513 legend. Please have a look at this as it appears it may have been an error when preparing the\n514 files for submission.\n516 visualize all 79 trees in one figure. Each line represents a tree, so this figure includes 79 trees.\n519 [Comment of Reviewer #2]\n521 a statistic to detect positively selected and introgressed genes (PSIGs)” Are these actually\n522 detecting genes or regions? You use both terms in the same sentence.\n524 statistic to identify frequently introgressed genes in the EAI genomes”\n525 [Comment of Reviewer #2]\n528 [Comment of Reviewer #2]\n530 genome assembly, it is NOT an accession for a gene set. I assume that you either used\n531 annotations from NCBI or Ensembl, both of which have many different versions which are\n532 regularly updated. You need to specify which annotations you used (NCBI or Ensembl) and\n533 exactly which version of the annotation.\n535 “As source of annotation, we used the source Bos taurus Annotation Release 106\n536 (GCF_002263795.1_ARS-UCD1.2_genomic.gtf) based on the NCBI assembly of\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have improved the manuscript considerably and have taken all my comments\ninto account.\nTherefore, I have no additional comments.\nReviewer #2 (Remarks to the Author):\nI have reviewed all of the material the authors provided. In my opinion, the authors have\naddressed the main issues that I have raised that needed to be addressed. At this point,\nthere seem to still be some issues regarding methodology, but those are unlikely to be\nresolved as they are differences of opinion on how something should be done. Likewise, I\nbelieve there may still be some differences of opinion regarding interpretation of results.\nHowever, I see nothing that would preclude this manuscript from being published, and any\nperceived issues addressed through the normal scientific process.\nRobert Schnabel"}]

Decision letter