Understanding the process and consequences of hybridization is one of the major challenges in evolutionary biology. A growing body of literature has reported evidence of ancient hybridization events or natural hybrid zones in primates, including humans; however, we still have relatively limited knowledge about the pattern and history of admixture because there have been little studies that simultaneously achieved genome-scale analysis and a geographically extensive sampling of wild populations. Our study applied double-digest restriction site-associated DNA sequencing to samples from the six localities in and around the provisional hybrid zone of rhesus and long-tailed macaques and evaluated population structure, phylogenetic relationships, demographic history, and geographic clines of morphology and allele frequencies. A latitudinal gradient of genetic components was observed, highlighting the transition from rhesus (north) to long-tailed macaque distribution (south) as well as the presence of one northern population of long-tailed macaques exhibiting unique genetic structure. Interspecific gene flow was estimated to have recently occurred after an isolation period, and the migration rate from rhesus to long-tailed macaques was slightly greater than in the opposite direction. Although some rhesus macaque-biased alleles have widely introgressed into long-tailed macaque-populations, the inflection points of allele frequencies have been observed as concentrated around the traditionally recognized interspecific boundary where morphology discontinuously changed; this pattern was more pronounced in the X-chromosome than in autosomes. Thus, due to geographic separation before secondary contact, reproductive isolation could have evolved, contributing to the maintenance of an interspecific boundary and species-specific morphological characteristics.

Data structure:
.
├── README.txt  # This file
│
│
├── fastsimcoal2 (compressed in zip format)  # The data for fastsimcoal2
│   ├── setting  # Estimation files (.est) and template files (.tpl)
│   │   ├── default  # Main analysis
│   │   ├── n_her  # Rhesus population size is fixed at 239704 (479408 genomes)
│   │   ├── n_xue  # Rhesus population size is fixed at 71000 (142000 genomes)
│   │   ├── p_2fold  # SFS projection size is two-hold (40, 50) 
│   │   └── p_half  # SFS projection size if half (10, 13)
│   │
│   └── sfs  # Multidimensional SFS provided by easySFS.py
│       ├── default_bootstrap  # Main analysis. 100 bootstrap samples and the original.
│       ├── p_2fold  # SFS projection size is two-hold (40, 50) 
│       ├── p_half  # SFS projection size if half (10, 13) 
│       ├── sub1  # Subset 1. Populations close to interspecific boundary (RH-BSS, RH-WTPMH, and LT-WHM) are excluded.
│       ├── sub2  # Subset 2. Populations close to and disproportionately-far-away from interspecific boundary (RH-BSS, RH-WTPMH, LT-WHM, and LT-Sumatra) are excluded.
│       └── sub3  # Subset 3. Captive populations (RH-China and LT-Sumatra) are excluded.
│
│
├── fineRADstructure (compressed in zip format)  # The haplotype data used in fineRADstructure (the oupput of Stacks populations command with --radpainter option)
│   ├── autosome.haps.radpainter  # All the samples
│   └── autosome_wild.haps.radpainter  # Wild samples (RH-China and LT-Sumatra are excluded)
│
│
├── geocline (compressed in zip format)  # The data and R script used in geographic cline analysis
│   ├── genome.csv  # Location information for SNPs data
│   ├── geocline_functions.R  # R script of supportive functions
│   ├── mtd.csv  # mtDNA data (after Figure 2 of Bunlungsup et al., 2017)
│   ├── rtl.csv  # Relative tail length data (after Table 1 of Hamada et al., 2015)
│   └── ych.csv  # Y-chromosome data (after Figure 2 of Bunlungsup et al., 2017)
│
│
├── nj (compressed in zip format)  # The input and output of PAUP
│   ├── nj_a.dist  # Uncorrected P-distance of autosomal SNPs calculated by PAUP. This is used for Neighbor-net analysis in SplitsTree4.
│   ├── nj_a.nexus  # Nexus of autosome used in PAUP NJ analysis
│   ├── nj_a.support.tre  # NJ tree with bootstrap support value of autosome
│   ├── nj_x.nexus  # Nexus of X-chromosome used in PAUP NJ analysis
│   ├── nj_x.support.tre  # NJ tree with bootstrap support value of X-chromosome
│   ├── nj_y.nexus  # Nexus of Y-chromosome used in PAUP NJ analysis
│   └── nj_y.support.tre  # NJ tree with bootstrap support value of Ychromosome
│
│
├── sex_pop_list.txt  # Tab-delimited text file of ID, sex, and population
│
│
└── vcf_filtered (compressed in zip format)  # Filtered vcf files
│   ├── autosome.vcf  # This is used in population structure analyses. This is also converted to nexus format for phylogenetic analyses using vcf2phylip.py.
│   ├── autosome_sfs.vcf  # This is used for making SFS (the input for easySFS.py).
│   ├── x_chromosome.vcf  # This is used in population structure analyses. This is also converted to nexus format for phylogenetic analyses using vcf2phylip.py.
│   └── y_chromosome.vcf  # This is converted to nexus format for phylogenetic analyses using vcf2phylip.py.
│
│
└── vcf_raw (compressed in zip format)  # Raw vcf files (the output of Stacks populations command with --vcf option)
    ├── autosome_raw.vcf
    ├── x_chromosome_raw.vcf
    └── y_chromosome_raw.vcf