Innovations in ancient DNA (aDNA) preparation and sequencing technologies have exponentially increased the quality and quantity of aDNA data extracted from ancient biological materials. The additional temporal component from the incoming aDNA data can provide improved power to address fundamental evolutionary questions like characterising selection processes that shape the phenotypes and genotypes of contemporary populations or species. However, utilising aDNA to study past selection processes still involves considerable hurdles such as how to eliminate the confounding effect of genetic interactions in the inference of selection. To circumvent this challenge, in this work we extend the method introduced by He et al. (2022) to infer temporally variable selection from the data on aDNA sequences with the flexibility of modelling linkage and epistasis. Our posterior computation is carried out through a robust adaptive version of the particle marginal Metropolis-Hastings algorithm with a coerced acceptance rate. Moreover, our extension inherits their desirable features like modelling sample uncertainties resulting from the damage and fragmentation of aDNA molecules and reconstructing underlying gamete frequency trajectories of the population. We assess the performance and show the utility of our procedure with an application to ancient horse samples genotyped at the loci encoding base coat colours and pinto coat patterns.