Recent advances in DNA sequencing technology, which reduced costs and increased throughput and accuracy, have driven subsequent advances in population genomics methods for detecting traces of natural selection in DNA fragments. In a recombining chromosome, a subgenomic region under natural selection typically exhibits different levels of variation and differentiation than the rest of the genome (Li & Ralph, 2019). Hence, it can be considered an anomaly that deviates from the overall population structure (François, Martins, Caye, & Schoville, 2016; Haasl & Payseur, 2016). Identifying such anomalies in molecular data is of great significance since the respective footprints of localized natural selection can provide insight into the adaptation process of a population to its environment through different generations. One of the most frequently used statistics to detect genomic regions under selection is the Fixation Index (F ST), which was introduced to quantify population differentiation based on the Wright-Fisher model (Wright, 1949). Several F ST variants are widely employed in population genomics (Hudson, Slatkin, & Maddison, 1992; Weir & Cockerham, 1984; Weir & Ott, 1997) because high F ST values can be an indication of local adaptation. However, when the population history deviates from the Wright-Fisher model, or when evolutionary history is described by a hierarchical population structure model, hypothesis testing becomes a challenge because the F ST distribution that accounts for the neutral demographic model of the population under study is not known. In this case, F ST-based methods that do