Your question is a bit vague, but let me try the following statement, which might be the kind of answer you are looking for: If $M$ is a manifold and $S\to M$ is a vector bundle over $M$ endowed with a connection $\nabla$, then a 'total differential equation' for sections $s$ of $S$ is an equation of the form $\nabla_X s = 0$ for all vector fields $X$. Then the Frobenius integrability conditions for this system are simply $\Omega^\nabla=0$, where $\Omega^\nabla$ is the curvature of $\nabla$, i.e., it is the $2$-form on $M$ with values in $\mathrm{End}(S)$ that is defined by
$$
\Omega^\nabla(X,Y)s = \nabla_X\nabla_Ys - \nabla_Y\nabla_Xs - \nabla_{[X,Y]}s
$$
for all vector fields $X$ and $Y$ on $M$. So, yes, the Frobenius integrability conditions are 'covariant', i.e., they can be formulated independently of coordinates.
Now, one may wonder, "But suppose I start with some linear total differential system
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{bi}(x) u^b \tag{1}
$$
in local coordinates. Does this apply?" The answer is 'yes': One regards the $u^a$ as the components of a section $s = (u^a)$ of a trivialized vector bundle $S$ and the matrices $\Gamma_i = \bigl(\Gamma^a_{bi}(x)\bigr)$ as sections of $\mathrm{End}(S)$ and one defines a connection $\nabla$ on this trivial bundle by the rule
$$
\nabla_{\frac{\partial}{\partial x^i}}s = \frac{\partial s}{\partial x^i} - \Gamma_i\ s.
$$
Then the curvature of $\nabla$ vanishes if and only if the Frobenius conditions are satisfied for the system (1).
If one has an inhomogeneous linear system
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{bi}(x) u^b + \Lambda^a_i(x), \tag{2}
$$
one can convert it to a homogeneous system by replacing $\Lambda^a_i(x)$ by $\Lambda^a_i(x)u^0$ in the above equation and adjoining the equations
$$
\frac{\partial u^0}{\partial x^i} = 0,
$$
making it a homogeneous linear system in one more variable. (Effectively, one allows the $a$ and $b$ indices (which run from, say, $1$ to $s$) to attain the value $0$, and sets $\Gamma^0_{bi}= \Gamma^0_{0i}=0$ and $\Gamma^a_{0i}=\Lambda^a_i$.) Again, the Frobenius conditions are equivalent to the vanishing of the curvature of the connection associated to the augmented system (now in $s{+}1$ unknowns rather than $s$ unknowns).
If one goes all the way to a fully nonlinear total differential system, say,
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{i}(x,u), \tag{3}
$$
then one can't (usually) reformulate this as the vanishing of the covariant derivatives of a section of something. However, there is still a 'covariant' coordinate-free interpretation of the Frobenius condition in terms of a tensor, but the tensor is defined on the total space $S$ (which has the $x$'s and $u$'s as local coordinates). The point is then, that, on $S$, one has a plane field $D$ (sometimes called a 'distribution') defined by the Pfaffian equations
$$
du^a - \Gamma^a_{i}(x,u)\ dx^i = 0
$$
(summation on $i$ intended). Associated to this plane field $D\subset TS$, there is a natural skew-symmetric linear operator
$$
\Phi: D\times D \to TS{\bigl/}D
$$
i.e., a section of the bundle $\bigl(TS{\bigl/}D\bigr)\otimes \Lambda^2(D^\ast)$ over $S$ that vanishes if and only if the system (3) satisfies the Frobenius conditions. It is this operator that is the 'covariant' interpretation of the Frobenius integrability conditions. Of course, this operator is defined via the Lie bracket. On sections $X$ and $Y$ of $D$, it is just
$$
\Phi(X,Y) \equiv [X,Y]\ \mathrm{mod}\ D,
$$
so one gets the usual statement that the Frobenius integrability conditions are equivalent to the condition that the sections of $D$ be closed under Lie bracket.
Finally, to get back to when one can interpret the Frobenius conditions as the vanishing of the curvature of a connection, there is a rather general formulation of this that goes back all the way to Lie: Suppose that the system (3) can be written in the form
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{i}(x,u)
= \gamma_{i\sigma}(x) F^{a\sigma}(u), \tag{4}
$$
for some functions $\gamma_{i\sigma}$ of $x$ and $F^{a\sigma}$ of $u$ (where summation over a new index $\sigma$ is assumed), and, further, assume that the vector fields
$$
U^\sigma = F^{a\sigma}(u)\ \frac{\partial\ \ }{\partial u^a}
$$
form a Lie algebra $L$, i.e., that there exist constants $c_\sigma^{\rho\tau}=-c_\sigma^{\tau\rho}$ such that $[U^\rho,U^\tau] = c_\sigma^{\rho\tau}U^\sigma$. Then the system (4) is what is known as a system of Lie type.
To avoid degenerate situations, one can assume, without loss of generality, that the $U^\sigma$ are a basis of $L$. In this case, one has that the Frobenius integrability condition for (4) is equivalent to the vanishing of the curvature of the $L$-valued connection $1$-form
$$
\omega = U^\sigma\ \gamma_{i\sigma}(x)\ dx^i.
$$
The use of this form of the equations belongs in the theory of symmetry analysis of differential equations, about which there is a very large literature.
Of course, recognizing when (3) can be written in the form (4) is not a trivial matter; there are methods for doing this involving Cartan's Method of Equivalence, but that would take a while to explain.
The answer to Q2 for $n=3$ is actually 'no, without some nondegeneracy hypotheses'. The reason is as follows:
The curvature tensor $\mathcal{R}= R_{ijkl}\,(\mathrm{d}x^i\wedge\mathrm{d}x^j)\circ(\mathrm{d}x^k\wedge\mathrm{d}x^l)$, with all its indices lowered, is a section of the subbundle $K(M)\subset S^2\bigl(\Lambda^2(T^*M)\bigr)$ that is the kernel of the natural linear mapping
$$
S^2\bigl(\Lambda^2(T^*M)\bigr)\longrightarrow \Lambda^4(T^*M).
$$
You are asking, for a given section $\mathcal{R}$, whether there exists a metric $g$ such that $\mathrm{Riem}(g) = \mathcal{R}$. There is no pointwise algebraic condition on the section $\mathcal{R}$ imposed by this equation (that was what Q1 was about), but, as you note, there is the second Bianchi identity
$$
\mu(\nabla^g\mathcal{R}) = 0,
$$
where $\mu:K(M)\otimes T^*M\to \Lambda^2(T^*M)\otimes\Lambda^3(T^*M)$ is the natural skewsymmetrization operation. Since $\nabla^g$ depends on one derivative of the metric $g$, the above equation with a given $\mathcal{R}$ can be regarded as a first-order system of equations on $g$. When $n=3$, this is at most $3$ equations, the rank of the bundle $\Lambda^2(T^*M)\otimes\Lambda^3(T^*M)$.
Now, for a point $p\in M$ satisfying $\mathcal{R}(p)=0$, the value $\nabla^g\mathcal{R}(p)$ does not depend on $g$ (just look at the formula for $\nabla^g\mathcal{R}$ in local coordinates). Thus, if $\mathcal{R}(p)=0$, but $\mu(\nabla^{g_0}\mathcal{R})(p) \not= 0$ for some metric $g_0$, then $\mu\bigl(\nabla^{g}\mathcal{R}(p)\bigr) \not= 0$ for all metrics $g$ and hence there is no open neighborhood of $p$ on which the equations $\mu\bigl(\nabla^{g}\mathcal{R}\bigr) = 0$ have a solution $g$. In particular, the original system $\mathrm{Riem}(g) = \mathcal{R}$ has no solution in a neighborhood of such a point $p$.
To convince yourself that such examples exist when $n=3$, just note that the rank of $K(M)=S^2\bigl(\Lambda^2(T^*M)\bigr)$ in this case is $6$, so an arbitrary section that vanishes at $p$ will have $6\times 3 = 18$ independent first derivatives. Thus, the map $\mu:K(M)\otimes T^*M\to \Lambda^2(T^*M)\otimes\Lambda^3(T^*M)$ is surjective (and the rank of the target bundle is $3$), so that the generic section of $K(M)$ that vanishes at $p$ will not satisfy the second Bianchi identity at $p$ for any metric $g$.
However, suppose that $n=3$ and that $\mathcal{R}$ is a nondegenerate section of $K(M)=S^2\bigl(\Lambda^2(T^*M)\bigr)$. I proved (back in the early 1980s) that, when $\mathcal{R}$ is real-analytic, there always exist local solutions to the equation $\mathrm{Riem}(g) = \mathcal{R}$. Specifically, I showed that, in this case, in addition to the $6$ second-order equations that these equations represent on $g$ and the $3$ first-order equations on $g$ that $\mu\bigl(\nabla^{g}\mathcal{R}\bigr) = 0$ represents, there is one more first-order equation $Q_\mathcal{R}(g)=0$ on $g$ that is satisfied by any metric $g$ that satisfies $\mathrm{Riem}(g) = \mathcal{R}$. Then I proved that the combined overdetermined system of $6$ second-order equations and $4$ first-order equations for $g$ is involutive, so that an application of the Cartan-Kähler Theorem proves local solvability. Unfortunately, the involutive system is never either hyperbolic or elliptic, though it can be of real principal type.
I never published my proof, but later, Dennis DeTurck and Deane Yang studied the overdetermined system that I wrote down and published a proof of its local solvability in the smooth category. See Deturck and Yang, Local existence of smooth metrics with prescribed curvature, Nonlinear problems in geometry (Mobile, Ala., 1985), 37–43,
Contemp. Math., 51, Amer. Math. Soc., Providence, RI, 1986.
Best Answer
Part of the difficulty in providing an answer to your question is the fact that the expression "the integrability condition" is a somewhat vague notion, and it's used in slightly different senses in different contexts.
The usual, somewhat imprecise, sense is that, for a given system of PDE, its 'integrability condition' is the set of necessary and sufficient conditions that it have local solutions. Slightly more precise usage would be, given a PDE system that depends on some data (in your case, the equations (1), where the data is the curvature $R$), to find the necessary and sufficient conditions on the data that ensure that the PDE system have local solutions.
In this sense, the Bianchi identities (your equation (2)) are usually not the integrability condition for the system (1) but are only part of the integrability condition, and they have to be interpreted appropriately. The point is that, when $R$ is specified, equation (2) is an algebraic equation $\Gamma$ must satisfy if it is to satisfy the first order PDE system (1). [It may not be immediately apparent that (2) is an algebraic equation on $\Gamma$, but if you look at the definition of $\nabla$, you'll see that the left hand side of (2) only involves the unknown coefficients of $\Gamma$ linearly with no differentiation of them, and the right hand side involves the unknown $S$, which is computed linearly from the unknown $\Gamma$.] Thus, a necessary condition on the data $R$ that (1) have a solution is that there be at least one solution $\overline\Gamma$ to the inhomogeneous linear algebraic system (2) (whose coefficients come from $R$ and its first derivatives). Generally, this necessary condition is not sufficient, though, as examples in dimensions $4$ and above show.
However, what Deane Yang's reference eudml.org/doc/74779 (DeTurck & Talvacchia, 1987) shows is that, when (i) $n=3$, (ii) $R$ is sufficiently generic (in an appropriate sense) while satisfying the integrability condition (2), and (iii) $R$ is real-analytic, then the system (1) is locally solvable. The system of PDE for $\Gamma$ that has to be solved, though, is highly nonlinear and, in a sense, overdetermined, so that the Cartan-Kähler Theorem or one of its modern versions has to be applied. There is no uniqueness and no way explicitly to solve the equations for $\Gamma$ for a given generic $R$ that satisfies the integrability condition. In that sense, your second 'reverse question' has no good answer.
Finally, let me remark that one could conceivably want to answer the more restrictive problem of specifying both the curvature $R$ and torsion $S$ of a connection $\Gamma$ on a manifold $M$. This problem is overdetermined even when $n=3$ and, generally, has no solution. This problem also has an additional set of integrability conditions of the form $\nabla S = F(R,S)$ (for an explicit expression $F(R,S)$ that one can write down easily) and which are also inhomogeneous linear algebraic in $\Gamma$, in addition to the Bianchi conditions $\nabla R = G(R,S)$ that you already have written down as (2). In general, these combined integrability conditions are not sufficient for the generic pair $(S,R)$ that satisfy them in order for a solution $\Gamma$ to exist.