Biochemical and evolutionary interactions between mitochondrial and nuclear genomes (‘mitonuclear interactions’) are proposed evolutionary drivers of sexual reproduction, sexual selection, adaptation, and speciation. We investigated the role of pre-mating isolation in maintaining functional mitonuclear interactions in wild populations bearing diverged proposed co-adapted mitonuclear genotypes. Two lineages of eastern yellow robin Eopsaltria australis—putatively climate-adapted to ’inland’ and ‘coastal’ climates—differ by ~7% of mitochondrial DNA positions, whereas nuclear genome differences are concentrated into a sex-linked region enriched with genes with mitochondrial functions. This pattern can be explained by female-linked selection accompanied by male-mediated gene flow across the narrow hybrid zone where the two lineages coexist. It remains unknown whether lineage divergence is driven by intrinsic incompatibilities (particularly in females, under Haldane’s rule), extrinsic incompatibilities, both, or other drivers. We tested whether lineage divergence could be facilitated by non-random mate-pairing with respect to partners’ mitolineage or nuclear Z sex-chromosome DNA sequences, which differ between the lineages. We used field-, Z-linked-, and mitolineage data from two locations where the lineages hybridize, to test whether females mate disproportionately with (1) males of their own mitolineage and/or bearing similar Z-linked variation, as might be expected if hybrids experience intrinsic incompatibilities, or (2) putatively locally-adapted males, as expected under environmental selection. Comparing field observations with simulations provided no evidence of non-random mating, thus the observed patterns consistent with reduced female gene flow likely occur post-mating. Future tests of female-biased mortality at different life stages and habitat selection may clarify any mechanisms of selection.

Identifying sex-linked markers in genomic datasets is important, because their analyses can reveal sex-specific biology, and their presence in supposedly neutral autosomal datasets can result in incorrect estimates of genetic diversity, population structure and parentage. But detecting sex-linked loci can be challenging, and available scripts neglect some categories of sex-linked variation. Here, we present new R functions to (1) identify and separate sex-linked loci in ZW and XY sex determination systems and (2) infer the genetic sex of individuals based on these loci. Two additional functions are presented, to (3) remove loci with artefactually high heterozygosity, and (4) produce input files for parentage analysis. We test these functions on genomic data for two sexually-monomorphic bird species, including one with a neo-sex chromosome system, by comparing biological inferences made before and after removing sex-linked loci using our function. We found that standard filters, such as low read depth and call rate, failed to remove up to 28.7% of sex-linked loci. This led to (i) overestimation of population FIS by ≤ 9%, and the number of private alleles by ≤ 8%; (ii) wrongly inferring significant sex-differences in heterozygosity, (iii) obscuring genetic population structure, and (iv) inferring ~11% fewer correct parentages. We discuss how failure to remove sex-linked markers can lead to incorrect biological inferences (e.g., sex-biased dispersal and cryptic population structure) and misleading management recommendations. For reduced-representation datasets with at least 15 known-sex individuals of each sex, our functions offer convenient, easy-to-use resources to avoid this, and to sex the remaining individuals.