If all conserved quantities of a system are known, can they be explained by symmetries?
There are already several good answers. However, the off-shell aspect related to Noether Theorem has not been addressed so far. (The words on-shell and off-shell refer to whether the equations of motion (e.o.m.) are satisfied or not.) Let us rephrase the problem as follows.
Consider a (not necessarily isolated) Hamiltonian system with $N$ degrees of freedom (d.o.f.). The phase space has $2N$ coordinates, which we denote $(z^1, \ldots, z^{2N})$.
(We shall have nothing to say about the corresponding Lagrangian problem.)
Symplectic structure. Usually, we work in Darboux coordinates $(q^1, \ldots, q^N; p_1, \ldots, p_N)$, with the canonical symplectic potential one-form
$$\vartheta=\sum_{i=1}^N p_i dq^i.$$ However, it turns out to be more efficient in later calculations, if we instead from the beginning consider general coordinates $(z^1, \ldots, z^{2N})$ and a general (globally defined) symplectic potential one-form $$\vartheta=\sum_{I=1}^{2N} \vartheta_I(z;t) dz^I,$$ with non-degenerate (=invertible) symplectic two-form $$\omega = \frac{1}{2}\sum_{I,J=1}^{2N} \omega_{IJ} \ dz^I \wedge dz^J = d\vartheta,\qquad\omega_{IJ} =\partial_{[I}\vartheta_{J]}=\partial_{I}\vartheta_{J}-\partial_{J}\vartheta_{I}. $$ The corresponding Poisson bracket is $$\{f,g\} = \sum_{I,J=1}^{2N} (\partial_I f) \omega^{IJ} (\partial_J g), \qquad \sum_{J=1}^{2N} \omega_{IJ}\omega^{JK}= \delta_I^K. $$Action. The Hamiltonian action $S$ reads $$ S[z]= \int dt\ L_H(z^1, \ldots, z^{2N};\dot{z}^1, \ldots, \dot{z}^{2N};t),$$ where $$ L_H(z;\dot{z};t)= \sum_{I=1}^{2N} \vartheta_I(z;t) \dot{z}^I- H(z;t) $$ is the Hamiltonian Lagrangian. By infinitesimal variation $$\delta S = \int dt\sum_{I=1}^{2N}\delta z^I \left( \sum_{J=1}^{2N}\omega_{IJ} \dot{z}^J-\partial_I H - \partial_0\vartheta_I\right)+ \int dt \frac{d}{dt}\sum_{I=1}^{2N}\vartheta_I \delta z^I, \qquad \partial_0 \equiv\frac{\partial }{\partial t},$$ of the action $S$, we find the Hamilton e.o.m. $$ \dot{z}^I \approx \sum_{J=1}^{2N}\omega^{IJ}\left(\partial_J H + \partial_0\vartheta_J\right) = \{z^I,H\} + \sum_{J=1}^{2N}\omega^{IJ}\partial_0\vartheta_J. $$ (We will use the $\approx$ sign to stress that an equation is an on-shell equation.)
Constants of motion. The solution $$z^I = Z^I(a^1, \ldots, a^{2N};t)$$ to the first-order Hamilton e.o.m. depends on $2N$ constants of integration $(a^1, \ldots, a^{2N})$. Assuming appropriate regularity conditions, it is in principle possible to invert locally this relation such that the constants of integration $$a^I=A^I(z^1, \ldots, z^{2N};t)$$ are expressed in terms of the $(z^1, \ldots, z^{2N})$ variables and time $t$. These functions $A^I$ are $2N$ (locally defined) constants of motion (c.o.m.), i.e., constant in time $\frac{dA^I}{dt}\approx0$. Any function $B(A^1, \ldots, A^{2N})$ of the $A$'s, but without explicit time dependence, will again be a c.o.m. In particular, we may express the initial values $(z^1_0, \ldots, z^{2N}_0)$ at time $t=0$ as functions $$Z^J_0(z;t)=Z^J(A^1(z;t), \ldots, A^{2N}(z;t); t=0)$$ of the $A$'s, so that $Z^J_0$ become c.o.m.
Now, let $$b^I=B^I(z^1, \ldots, z^{2N};t)$$ be $2N$ independent (locally defined) c.o.m., which we have argued above must exist. OP's title question in this formulation then becomes if there exist $2N$ off-shell symmetries of the (locally defined) action $S$, such that the corresponding Noether currents are on-shell c.o.m.?
Remark. It should be stressed that an on-shell symmetry is a vacuous notion, because if we vary the action $\delta S$ and apply e.o.m., then $\delta S\approx 0$ vanishes by definition (modulo boundary terms), independent of what the variation $\delta$ consists of. For this reason we often just shorten off-shell symmetry into symmetry. On the other hand, when speaking of c.o.m., we always assume e.o.m.
Change of coordinates. Since the action $S$ is invariant under change of coordinates, we may simply change coordinates $z\to b = B(z;t)$ to the $2N$ c.o.m., and use the $b$'s as coordinates (which we will just call $z$ from now on). Then the e.o.m. in these coordinates are just $$\frac{dz^I}{dt}\approx0,$$ so we conclude that in these coordinates, we have $$ \partial_J H + \partial_0 \vartheta_J=0$$ as an off-shell equation. [An aside: This implies that the symplectic matrix $\omega_{IJ}$ does not depend explicitly on time, $$\partial_0\omega_{IJ} =\partial_0\partial_{[I}\vartheta_{J]}=\partial_{[I} \partial_0\vartheta_{J]}=-\partial_{[I}\partial_{J]} H=0.$$ Hence the Poisson matrix $\{z^I,z^J\}=\omega^{IJ}$ does not depend explicitly on time. By Darboux Theorem, we may locally find Darboux coordinates $(q^1, \ldots, q^N; p_1, \ldots, p_N)$, which are also c.o.m.]
Variation. We now perform an infinitesimal variation $\delta= \varepsilon\{z^{I_0}, \cdot \}$, $$\delta z^J = \varepsilon\{z^{I_0}, z^J\}=\varepsilon \omega^{I_0 J},$$ with Hamiltonian generator $z^{I_0}$, where $I_0\in\{1, \ldots, 2N\}$. It is straightforward to check that the infinitesimal variation $\delta= \varepsilon\{z^{I_0}, \cdot \}$ is an off-shell symmetry of the action (modulo boundary terms) $$\delta S = \varepsilon\int dt \frac{d f^0}{dt}, $$ where $$f^0 = z^{I_0}+ \sum_{J=1}^{2N}\omega^{I_0J}\vartheta_J.$$ The bare Noether current is $$j^0 = \sum_{J=1}^{2N}\frac{\partial L_H}{\partial \dot{z}^J} \omega^{I_0 J}=\sum_{J=1}^{2N}\omega^{I_0J}\vartheta_J,$$ so that the full Noether current $$ J^0=j^0-f^0=-z^{I_0} $$ becomes just (minus) the Hamiltonian generator $z^{I_0}$, which is conserved on-shell $\frac{dJ^0}{dt}\approx 0$ by definition.
So the answer to OP's title question is Yes in the Hamiltonian case.
See also e.g. this, this & this related Phys.SE posts.
Yes, this is the opposite of Noether's theorem. So let's call our conserved quantity $A$ (we will consider just one conserved quantity for starters) and begin with $\left \{H, A \right \} = 0$ law for conservation. Because of the connection between Poisson bracket with flows on the phase space this tells you both that $\mathcal{L}_{V_H} A$ = 0 ($A$ is conserved in time evolution) and $\mathcal{L}_{V_A} H = 0$ (Hamiltonian is conserved under symmetry that has fundamental field $V_A$) as long as the vectors fields associated to the functions on phase space as $V_X = ({\rm d} X)^{\sharp}$ are non-degenerate in certain sense. Note that raising $(\cdot)^{\sharp}$ operator is defined obviously using the inverse symplectic form. This e.g. means that for example genuine constants (which are certainly also constants of motion) won't work because ${\rm d} C = 0$ and we obtain zero vector field.
On the other hand, as long as all the conserved quantities are non-degenerate we can always find the associated symmetry flows by integration of the above mentioned vector fields. But note that what we receive in the end are symmetries of the phase space. Whether these are also directly linked to some symmetry of the underlying position space (if there is such a space, that is; in general there need not be one but in usual applications we take the phase space of the position manifold $M$ as the cotangent bundle $T^*M$) is a question for further investigation. I'll try looking into it later, if I find some time.
I'll start with some terminological subtleties, that have to be accounted for, when it comes to "conserved quantities" or "integrals of motion".
First of all it is important to state on which variables the quantities can depend. In the field of differential equations, Hamiltonian dynamics and dynamical systems, a conserved quantity by definition does not explicitly depends on time -- it is usually a function, defined on a space (a manifold), that fully describes the state of your system (uses all the degrees of freedom as you have defined them).
Another thing to be taken into account -- It seems to be possible to take a function of a constant of motion, therefore arriving at another quantity that is constant, therefore producing any amount of constants of motion. So there is an implicit condition for functional independence between those $N$ constants of motion. Which can be formulated as the vanishing of the Poisson brackets between any pair of the quantities. Here is the reference about it on Wikipedia, also I've asked a question about it some time ago...
So I think that your implicit statement:
If a system has N degrees of freedom (DOF) and therefore N conserved independent quantities
Is not totally correct, if one takes into account what was said about the notion of "conserved quantity".
Concerning your question about finding the symmetry transformation, that corresponds to the quantity -- I've addressed the issue in this question.
Briefly -- given a "generator" $\delta G$ and some quantity $A$. The small transformation, generated by $\delta G$ is: $$A \to A+\delta A,\quad\quad \delta A = -\{\delta G, A\}$$
Putting $A = H$, and noting that, if $\delta G$ is a constant of motion, then $\{\delta G, H\} = 0$. One immediately arrives at the conclusion that the transformation, generated by $\delta G$ is the symmetry transformation. (Some examples are here.)