infer skyline plots from this type of data (e.g., Allen et al., 2012 ; Molfetti et al., 2013 ; Minhós et al., 2016 ). Also for microsatellite data, a composite-likelihood approach has been developed  (R package VarEff, Nikolic & Chevalet, 2014 ). Similar piecewise models to infer historical population sizes through time have been proposed in the context of population genomics (e.g., Li & Durbin, 2011 ; Terhorst, Kamm & Song, 2016 ). The methods discussed above assume a set of independent (unlinked) genetic markers. However, if a large proportion of the genome has been sequenced, the studied polymorphism are in linkage disequilibrium. Methods such as the Pairwise Sequentially Markovian Coalescent (PSMC, Li & Durbin, 2011 ) and its successors profit from the additional information of linkage disequilibrium for the inference. We will not further discuss this family of methods, as our focus here is on datasets of independent molecular markers, such as microsatellites, which remain reliable markers for low-budget projects. Note, however, the PSMC-like implementation on ABC by  Boitard et al. (2016) . The use of the skyline plot in the ABC framework was first proposed in  Burgarella et al. (2012) . Here, we provide a suite of R scripts (DIYABCskylineplot) to produce approximate-Bayesian-computation skyline plots from microsatellite data and evaluate its performance on simulated pseudo-data. We show the method to be useful for detecting population decline and expansion and discuss its limits. ABC skyline plots are then built for four study cases (whale shark, leatherback turtle, Western black-and-white colobus and Temminck’s red colobus) and compared with the demographic inference obtained by an alternative full likelihood method. Methods ABC skyline plot For a demographic skyline plot analysis within the ABC framework, our model consisted of a single population with constant size that instantaneously changes to a new size n times through time. The parameters (from present