In this article, we study the regularity of solutions to inhomogeneous time-fractional evolution equations involving anisotropic non-local operators in mixed-norm Sobolev spaces of variable order, with non-trivial initial conditions. The primary focus is on space-time non-local equations where the spatial operator is the infinitesimal generator of a vector of independent subordinate Brownian motions, making it the sum of subdimensional non-local operators. A representative example of such an operator is $(Δ_{x})^{β_{1}/2}+(Δ_{y})^{β_{2}/2}$. We establish existence, uniqueness, and precise estimates for solutions in corresponding Sobolev spaces. Due to singularities arising in the Fourier transforms of our operators, traditional methods involving Fourier analysis are not directly applicable. Instead, we employ a probabilistic approach to derive solution estimates. Additionally, we identify the optimal initial data space using generalized real interpolation theory.