Vertex-domain and temporal-domain smoothness of time-varying graph signals are cardinal properties that can be exploited for effective graph signal reconstruction from limited samples. However, existing approaches are not directly applicable when the signal's frequency occupancy changes with time. Moreover, while e.g., sensor network applications can benefit from directed graph models, the non-orthogonality of the graph eigenvectors can challenge spectral-based signal reconstruction algorithms. In this context, here we consider $K$-sparse time-varying signals with unknown frequency supports. By exploiting the smoothness of the varying graph frequency supports and employing shift operations over directed graphs, we study joint sampling of multiple varying signals based on Schur decomposition to reconstruct each signal by orthogonal frequency components. Firstly, joint frequency support of the multiple signals is identified by proposing a two-stage Individual-Joint sampling scheme. Based on the estimated frequency support, the GFT coefficients of each signal can be recovered using data collected in individual sampling stage. Greedy algorithms are proposed for vertex set selection and graph shift order selection, which enable a robust signal reconstruction against additive noise. Considering the signals in applications may be approximately $K$-sparse, we further exploit the samples in both individual and joint sampling stages and investigate the optimal signal reconstruction as a convex optimization problem with adaptive frequency support selection. The proposed optimal sampling and reconstruction algorithms outperform several existing schemes in random network and sensor network data gathering.