We study the problem of characterizing and computing the Gaussian nonanticipative rate-distortion function (NRDF) of partially observable multivariate Gauss-Markov processes with mean-squared error (MSE) distortion constraints. First, we extend Witsenhausen's "tensorization" approach originally used for single-letter random variables to causal processes, to obtain a new modified representation of the NRDF for the specific problem. For time-varying vector processes, we prove conditions so that the new modified NRDF is achieved and study its implications when it is not achievable. For both cases (which correspond to different bounds), we derive the characterization and the optimal realization, whereas we give the optimal numerical solution using the semidefinite programming (SDP) algorithm. Interestingly, the realization (for both bounds) is shown to be a linear functional of the current time-sufficient statistic of the past and current observations signals. For the infinite time horizon, we give conditions to ensure the existence of a time-invariant characterization from the finite-time horizon problems and a numerical solution using the SDP algorithm. For the time-invariant characterization, we also give strong structural properties that enable an optimal and approximate solution via a reverse-waterfilling algorithm implemented via an iterative scheme, which executes much faster than the SDP algorithm. For both finite and infinite time horizons, we study the special case of scalar processes. Our results are corroborated with various simulation studies and are also compared with existing results in the literature.