Let's address the questions in turn.
The first one concerns the bivariate Normal distribution of $(X_1,X_2).$ This is determined by its mean $(\mu_1,\mu_2),$ the variances $\sigma_1\gt 0$ and $\sigma_2\gt 0,$ and the correlation coefficient $\rho = \rho_{12}.$ The theory of least-squares regression in the bivariate Normal setting teaches us that the distribution of $Z_1=(X_1-\mu_1)/\sigma_1$ conditional on $Z_2=(X_2-\mu_2)/\sigma_2$ is Normal with mean $\rho Z_2.$ Since $Z_2$ has a standard Normal distribution, we may directly compute
$$\begin{aligned}
E[X_1\mid a \le X_2 \le b] &= E[\mu_1+\sigma_1 Z_1\mid a \le \mu_2 + \sigma_2 Z_2 \le b]\\
&= \mu_1 + \sigma_1\,E[Z_1 \mid (a-\mu_2)/\sigma_2 \le Z_2 \le (b-\mu_2)/\sigma_2]\\
&= \mu_1 + \sigma_1 \,E[\rho Z_2 \mid (a-\mu_2)/\sigma_2 \le Z_2 \le (b-\mu_2)/\sigma_2]\\
&= \mu_1 + \sigma_1\rho \,E[Z_2 \mid (a-\mu_2)/\sigma_2 \le Z_2 \le (b-\mu_2)/\sigma_2].
\end{aligned}$$
Abbreviate things by setting
$$\alpha=(a-\mu_2)/\sigma_2\text{ and }\beta=(b-\mu_2)/\sigma_2$$
for the standardized interval endpoints. Letting $\phi$ be the standard Normal PDF and recalling that $\mathrm{d}\phi(z) = -z\,\phi(z)\mathrm{d}z,$ it is easy to compute this conditional expectation as
$$\begin{aligned}
E[Z_2 \mid \alpha\le Z_2 \le \beta] &= \frac{1}{\Phi(\beta)-\Phi(\alpha)}\int_\alpha^\beta z \phi(z)\,\mathrm{d}z \\
&= \frac{1}{\Phi(\beta)-\Phi(\alpha)}\int_\alpha^\beta -\mathrm{d}\phi(z) \\
&= \frac{\phi(\alpha)-\phi(\beta)}{\Phi(\beta)-\Phi(\alpha)}.
\end{aligned}$$
This result suggests we should tackle the trivariate problem by first making similar simplifications: that is, standardize the variables and perform linear regression. Generalize the preceding notation to include a third subscript for the third variable. As before,
$$\begin{aligned}
E[X_1\mid a_2 \le X_2 \le b_2,\, a_3 \le X_3 \le b_3] &= E[\mu_1+\sigma_1 Z_1\mid a_2 \le \mu_2 + \sigma_2 Z_2 \le b_2, \ldots]\\
&= \mu_1 + \sigma_1\,E[Z_1 \mid (a_2-\mu_2)/\sigma_2 \le Z_2 \le (b_2-\mu_2)/\sigma_2, \ldots].
\end{aligned}$$
(This handles the case where the value of $(X_2,X_3)$ is conditioned on a rectangle rather than a square; it's no more difficult to solve and the notation makes the pattern a little clearer.)
The regression of $Z_1$ on $(Z_2,Z_3)$ is found by finding coefficients $(\gamma_2,\gamma_3)$ that minimize $$E[(Z_1 - \gamma_2Z_2 - \gamma_3Z_3)^2] = 1 - 2\rho_{12}\gamma_2 - 2\rho_{13}\gamma_3 + 2\rho_{23}\gamma_2\gamma_3 + \gamma_2^2 + \gamma_3^2.$$ Assuming $(X_1,X_2,X_3)$ has a nondegenerate distribution, its unique critical point (where the gradient vanishes), which must be the global minimum, therefore occurs when
$$\pmatrix{\rho_{12}\\\rho_{13}} = \pmatrix{1 & \rho_{23} \\ \rho_{23} & 1}\pmatrix{\gamma_2\\\gamma_3}$$
with solution
$$\pmatrix{\gamma_2\\\gamma_3} = \frac{1}{1-\rho_{23}^2} \pmatrix{\rho_{12} - \rho_{23} \rho_{13} \\ \rho_{13} - \rho_{23}\rho_{12}}.$$
Thus $$E[Z_1 \mid (Z_2,Z_3)] = \gamma_2 E[Z_2 \mid (Z_2,Z_3)]+ \gamma_3 E[Z_3\mid (Z_2,Z_3)]$$ and we may proceed as before to find this via integration. Again, to simplify the notation, write
$$\alpha_i = (a_i - \mu_i)/\sigma_i,\ \beta_i = (b_i - \mu_i)/\sigma_i$$
for $i=2, 3,$ and define the standard bivariate Normal conditional expectation functions
$$\begin{aligned}
\Psi_i([a,b],\,[c,d]\mid \rho) &= E[Z_i\mid Z_2\in[a,b],\ Z_3\in[c,d]] \\
&= \frac{\iint_{[a,b]\times[c,d]} z_i\,\phi(z_2,z_3\mid \rho)\,\mathrm{d}z_2\mathrm{d}z_3}{\iint_{[a,b]\times[c,d]} \phi(z_2,z_3\mid \rho)\,\mathrm{d}z_2\mathrm{d}z_3}
\end{aligned}$$
where
$$\phi(z_2,z_3\mid \rho) = \frac{1}{2\pi(1-\rho^2)} \exp\left(-\frac{z_2^2-2\rho z_2z_3+z_3^2}{2(1-\rho^2)}\right)$$
is the standard bivariate Normal PDF (for variables with correlation $\rho$). Our result (at the top of this section) then is
$$\begin{aligned}
&E[X_1\mid a_2 \le X_2 \le b_2,\, a_3 \le X_3 \le b_3] \\
&= \mu_1 + \sigma_1\,E[Z_1 \mid (a_2-\mu_2)/\sigma_2 \le Z_2 \le (b_2-\mu_2)/\sigma_2, \ldots] \\
&= \mu_1 + \sigma_1 \left(\gamma_2 \Psi_2([\alpha_2,\beta_2],[\alpha_3,\beta_3]\mid\rho_{23}) + \gamma_3 \Psi_3([\alpha_2,\beta_2],[\alpha_3,\beta_3]\mid\rho_{23})\right).
\end{aligned}$$
The integrals in the definition of the $\Psi_i$ are over a rectangle within the domain of the bivariate Normal distribution of $(Z_2,Z_3).$ Unless $\rho_{23}=0$ (when $Z_2$ and $Z_3$ are independent) the formula for this is messy; unless you really need a formula for further analysis, numerical integration may be most appropriate.
