Given $f:\mathbb R^n\to \mathbb R^m$, let $Df:\mathbb R^n\to \mathbb R^{m\times n}$ denote its Jacobian. Writing $f$ in terms of its $m$ component functions $f_1,f_2,\dots,f_m:\mathbb R^n\to \mathbb R$, we see that
$$
D\begin{bmatrix}f_1\\\vdots\\f_m\end{bmatrix}=\begin{bmatrix}Df_1\\\vdots\\Df_m\end{bmatrix}
$$
i.e. the rows of $Df$ are the Jacobians of the components of $f$. Therefore, in order to solve the equation $Df=F$, it suffices to solve each of the equations $Df_k=F_k$ seperately, for each $k\in \{1,\dots,m\}$, where $F_k$ is the $k^{th}$ row of $F$. In other words, we can restrict our attention to the $m=1$ case.
In this case, the Jacobian is just the gradient, so the question becomes
Given $F:\mathbb R^n\to \mathbb R^n$, when does there exist $f:\mathbb R^n\to \mathbb R$ for which $\nabla f=F$?
That is, how can we tell when a vector field is conservative? As long as $F$ is differentiable, with continuous partial derivatives, an obvious necessary condition is
$$
\forall i,j\in \{1,\dots,n\}:\frac{\partial F_j}{\partial x_i}=\frac{\partial F_i }{\partial x_j}
$$
The necessity of this condition follows from Schwarz's theorem, since if $F=\nabla f$, then
$$
\frac{\partial F_j}{\partial x_i}
=\frac{\partial }{\partial x_i}\frac{\partial f}{\partial x_j}
=\frac{\partial }{\partial x_j}\frac{\partial f}{\partial x_i}
=\frac{\partial F_i}{\partial x_j}.
$$
Schwarz's theorem requires $f$ to have continuous second partial derivatives, which is why I included the differentiability condition on $F$.
It turns out this condition is sufficient as well. In fact, it remains true for function $F:E\to \mathbb R^n$, where $E\subseteq \mathbb R^n$ is open, as long as $E$ is simply connected, meaning it does not have any "holes." However, I cannot prove this fact here, since it requires developing de Rham cohomology.
The maximum of any finite number of continuous functions is also continuous, cf. proof.
You correctly identified a symmetry for evaluating this integral that can be formalized as follows: For each permutation $\tau\in S_n$ (here, $S_n$ denotes the group of permutations on $\{1,\dots,n\}$), consider the Borel set $A_\tau = \{(x_1,\dots,x_n)\in[0,1]^n: x_{\tau(1)}\ge x_{\tau(2)}\ge\dots\ge x_{\tau(n)}\}$. We have $$\bigcup_{\tau\in S_n} A_\tau = [0,1]^n$$ and the pairwise intersection of any two distinct $A_\tau$ is a set of Lebesgue measure $0$ (exercise).
Therefore,
$$\int_{[0,1]^n} f(x)\,\mathrm dx = \sum_{\tau\in S_n}\int_{A_n} f(x)\,\mathrm dx.$$
Now notice, using the change of variables $(x_1,\dots,x_n)\mapsto(x_{\tau(1)},\dots, x_{\tau(n)})$ (details left to you), that for all $\tau\in S_n$, $$\int_{A_\tau} f(x)\,\mathrm dx = \int_{A_{\text{id}}} f(x)\,\mathrm dx,$$ where $\text{id}\in S_n$ denotes the identity map, i.e. the permutation leaving all elements fixed.
$S_n$ contains $n!$ elements and $\int_{A_{\text{id}}} f$ is exactly the integral you wrote down with slightly different notation.
The result $\frac{n-1}{n+1}$ is also correct. Alternatively you could note that if $X_1,\dots, X_n$ are independent, $\text{Uniform}([0,1])$-distributed random variables, then $$\mathbb E\left(\sup_{k\in\{1,\dots,n\}} X_k\right) = \frac{n}{n+1}$$ and $$\mathbb E\left(\inf_{k\in\{1,\dots,n\}} X_k\right) = \frac{1}{n+1}$$ so $$\mathbb E\left(\sup_{k\in\{1,\dots,n\}} X_k-\inf_{k\in\{1,\dots,n\}} X_k\right) = \frac{n-1}{n+1}.$$
Best Answer
We can rewrite the integral as
$$\int\limits_{[0,1]^n\cap\sum_i x_i \leq 1}\sum_i x_i \ d^n \mathbf{x} + \int\limits_{[0,1]^n\cap\sum_i x_i \geq 1}d^n \mathbf{x} = 1 - \int\limits_{[0,1]^n\cap\sum_i x_i \leq 1}1-\sum_i x_i \ d^n \mathbf{x}$$
$$ = 1 - \int\limits_{[0,1]^{n+1}\cap\sum_i x_i \leq 1} d^{n+1} \mathbf{x} \equiv 1 - \operatorname{vol}(n+1)$$
in other words, $1$ minus the volume of the sliced corner in $n+1$ dimensions. You can show that
$$\operatorname{vol}(n+1) = \int\limits_0^1\cdots \int\limits_0^{1-x_1-\cdots-x_{n-1}}(1-x_1-\cdots-x_n) \:dx_n\cdots dx_1$$ $$ = \int\limits_0^1\cdots\int\limits_0^{1-x_1-\cdots-x_{n-2}}\frac{(1-x_1-\cdots-x_{n-1})^2}{2}\:dx_{n-1}\cdots dx_1$$
$$= \cdots = \int_0^1 \frac{(1-x_1)^n}{n!}\:dx_1 = \frac{1}{(n+1)!}$$
Therefore the integral in question is $\boxed{1-\frac{1}{(n+1)!}}$