In this paper we study the system $$\begin{aligned}&\min \biggl \{-\mathcal H u_i(x,t)-\psi _i(x,t),u_i(x,t)-\max _{j\ne i}(-c_{i,j}(x,t)+u_j(x,t))\biggr \}=0,\\&u_i(x,T)=g_i(x),\ i\in \{1,\ldots ,d\}, \end{aligned}$$ where \((x,t)\in \mathbb R ^{N}\times [0,T]\). A special case of this type of system of variational inequalities with terminal data occurs in the context of optimal switching problems. We establish a general comparison principle for viscosity sub- and supersolutions to the system under mild regularity, growth, and structural assumptions on the data, i.e., on the operator \(\mathcal H \) and on continuous functions \(\psi _i\), \(c_{i,j}\), and \(g_i\). A key aspect is that we make no sign assumption on the switching costs \(\{c_{i,j}\}\) and that \(c_{i,j}\) is allowed to depend on \(x\) as well as \(t\). Using the comparison principle, the existence of a unique viscosity solution \((u_1,\ldots ,u_d)\) to the system is constructed as the limit of an increasing sequence of solutions to associated obstacle problems. Having settled the existence and uniqueness, we subsequently focus on regularity of \((u_1,\ldots ,u_d)\) beyond continuity. In this context, in particular, we assume that \(\mathcal H \) belongs to a class of second-order differential operators of Kolmogorov type of the form: $$\begin{aligned} \mathcal H =\sum _{i,j=1}^m a_{i,j}(x,t)\partial _{x_i x_j}+\sum _{i=1}^m a_i(x,t)\partial _{x_i} +\sum _{i,j=1}^N b_{i,j}x_i\partial _{x_j}+\partial _t, \end{aligned}$$ where \(1\le m\le N\). The matrix \(\{a_{i,j}(x,t)\}_{i,j=1,\ldots ,m}\) is assumed to be symmetric and uniformly positive definite in \(\mathbb R ^m\). In particular, uniform ellipticity is only assumed in the first \(m\) coordinate directions, and hence, \(\mathcal H \) may be degenerate.