In multivariate spatio-temporal analysis, we are faced with the formidable challenge of specifying a valid spatio-temporal cross-covariance function, either directly or through the construction of processes. This task is difficult as these functions should yield positive definite covariance matrices. In recent years, we have seen a flourishing of methods and theories on constructing spatio-temporal cross-covariance functions satisfying the positive definiteness requirement. A subset of those techniques produced spatio-temporal cross-covariance functions possessing the additional feature of nonstationarity. Here we provide a review of the state-of-the-art methods and technical progress regarding model construction. In addition, we introduce a rich class of multivariate spatio-temporal asymmetric nonstationary models stemming from the Lagrangian framework. We demonstrate the capabilities of the proposed models on a bivariate reanalysis climate model output dataset previously analyzed using purely spatial models. Furthermore, we carry out a cross-validation study to examine the advantages of using spatio-temporal models over purely spatial models. Finally, we outline future research directions and open problems.