Does the conservation of the Wronskian follow from Noether's principle? Noether's principle is the paradigm that symmetries of Hamiltonian and Lagrangian systems correspond to conservation laws of various kinds. Consider a one-dimensional harmonic oscillator
$$\tag{*} \ddot{x} + x =0.$$
If $x_1, x_2$ are two solutions, we have the conservation of their Wronskian:
$$
\frac d {dt}\left( x_1\dot{x}_2 -x_2\dot{x}_1\right) =0.$$

Does this conservation law correspond to a symmetry of the system? If yes, which one? 

A more general question is the following. Consider a Hamiltonian system, that is, the following ODE 
$$ \tag{**}
\frac{ d}{dt} \begin{bmatrix} q \\ p \end{bmatrix} = J\nabla_{q,p} H(q, p; t), $$ 
where $J$ is an antisymmetric matrix. Here one has Liouville's theorem which states that, if $\Omega(t)$ is a region of phase space that evolves with the flow of $(**)$, then 
$$\frac d{dt}\text{Vol }\Omega(t) = \frac d{dt} \iint_{\Omega(t)}dqdp=0.$$ (The aforementioned conservation of the Wronskian is a special case of this theorem, obtained by taking $\Omega(t)$ equal to the parallelogram spanned by $(x_1, \dot{x}_1)$ and $(x_2, \dot{x}_2)$).  
Same question as before: 

does this conservation law follow from a symmetry, and if yes, which one? 

The linked Wikipedia's page suggests that the conservation of volume in phase space follows from time translation invariance. This does not seem to me to be the case, because Liouville's theorem holds even in case of time-dependent Hamiltonians. The simplest example is a time dependent harmonic oscillator $\ddot{x}+b(t) x =0$. Here one still has conservation of the Wronskian:
$$\frac d{dt}\left(x_1\dot x_2 - x_2 \dot x_1\right) =-b(t)x_1x_2 +b(t)x_1x_2 =0.$$
 A: *

*First let us discuss the conserved
Wronskian
$$W(q_1,q_2)~=~q_1\dot{q}_2-q_2\dot{q}_1\tag{1}$$
for 1D systems. Within the 1D case, the property of a conserved Wronskian (1) doesn't seem to generally hold beyond a harmonic oscillator with explicit time dependence, i.e. Hooke's law
$$m\ddot{q}~\approx~-k(t)q ,\tag{2}$$
where the spring constant $k(t)$ may depend explicitly on time. (The $\approx$ symbol means equality modulo eom.)
In that case it is tempting to view the two solutions $q_1$ and $q_2$ as occurring along two perpendicular axes (which we will call $q_1$ and $q_2$, respectively, for simplicity) in the plane $\mathbb{R}^2$. (The plane $\mathbb{R}^2$ may be identified with the complex plane $\mathbb{C}$, cf. Valter Moretti's answer.) In other words, we are considering the 2D Lagrangian
$$ L_2~:=~\frac{m}{2}(\dot{q}_1^2+\dot{q}_2^2)-\frac{k(t)}{2}(q_1^2+q_2^2).\tag{3} $$
This Lagrangian (3) has a rotational symmetry, which by Noether's theorem means that the angular momentum
$$ L_3~:=~m(q_1\dot{q}_2-q_2\dot{q}_1)~=~m W(q_1,q_2),\tag{4}$$
and thereby the Wronskian (1) is conserved in time.


*Next, we should mention that OP's starting point is closely related to the covariant Hamiltonian formulation, cf. e.g. Ref. 1 and this & this Phys.SE posts. See also this Phys.SE post, which also starts with a Wronskian-like construction. (Below we are going to use Grassmann-odd variables, but it can equivalently be rephrased in the language of exterior calculus and wedge-products.) We start from the action
$$ S_0[q]~=~\int\! dt~L_0,\qquad L_0~:=~\frac{m}{2}\dot{q}^2- V(q),\tag{5} $$
with Euler-Lagrange (EL) equation
$$ 0~\approx~\frac{\delta S_0}{\delta q}~=~-m\ddot{q}-V^{\prime}(q). \tag{6}$$
We can introduce a Grassmann-odd nilpotent transformation$^1$
$${\rm s} q(t) ~=~ c(t), \qquad{\rm s}^2~=~0.\tag{7}$$
Consider now the Grassmann-odd action
$$\begin{align} S_1[q,c]~:=~&{\rm s}S_0~=~\int\! dt~L_1,\cr 
L_1~:=~&\frac{d}{dt}\left(m\dot{q}c\right)+ \frac{\delta S_0}{\delta q} c
~=~ m\dot{q}\dot{c}- V^{\prime}(q)c,\end{align}\tag{8} $$
with EL equations
$$\begin{align} 0~\approx~&\frac{\delta S_1}{\delta q}~=~-m\ddot{c}-V^{\prime\prime}(q)c, \cr 
0~\approx~&\frac{\delta S_1}{\delta c}~=~-m\ddot{q}-V^{\prime}(q). \end{align}\tag{9}$$
(The total time derivative term in eq. (8) is introduced for technical reasons to avoid higher-order time-derivatives in the Lagrangian $L_1$.) The Grassmann-odd action $S_1$ possesses by construction the Grassmann-odd symmetry (7), because ${\rm s}$ is nilpotent, i.e. squares to zero. We can therefore use a super-version of Noether's theorem to conclude that the corresponding Noether charge
$$ Q~=~c\frac{\partial L_1}{\delta \dot{q}} ~=~m c\dot{c} \tag{10}$$
is conserved on-shell
$$ \frac{dQ}{dt}~\approx~0. \tag{11} $$
The Noether charge (10) is the point-mechanical version of the symplectic 2-form current in Ref. 1. It constitutes the first step in a covariant Legendre transformation from Lagrangian to Hamiltonian formulation. (The word covariance refers to the fact that time and space are treated on equal footing. OP is only considering point mechanics, where covariance is not visible, but manifest Lorentz covariance becomes an issue in field theory. We leave it to the reader to generalize the above construction to field theory.)


*Finally, let us discuss Liouville's theorem. There are several versions of Liouville's theorem.
One version states that a Hamiltonian vector field $X_H=\{H,\cdot\}_{PB}$ on a symplectic manifold $(M,\omega)$ is divergence-free
$$ 0~=~{\rm div} X_{H}~=~\sum_{i=1}^{2n} \frac{\partial X_H^i}{\partial
z^i},\tag{12}$$ where $z^i$ are Darboux
coordinates/canonical
coordinates. Eq. (12) is hardwired into symplectic geometry. It follows from the fact that Hamiltonian vector fields preserve the symplectic structure
$$ {\cal L}_{X_H}\omega~=~ 0.\tag{13}$$
Hamilton's equations read
$$ \dot{z}^i~\approx~\{z^i,H\}_{PB}~=~-X_H^i.\tag{14}$$
If we form the current
$$V^{\mu}~=~(V^0,V^i)~=~(1,-X_H^i)~\approx~(\dot{z}^0,\dot{z}^i)~=~\dot{z}^{
\mu},\tag{15}$$
we can alternatively write (12) as
$$\sum_{\mu=0}^{2n}\frac{\partial V^{\mu}}{\partial z^{\mu}} ~=~0.\tag{16} $$
Here we have introduced the notation $z^0\equiv t$ for time.
Another version of Liouville's theorem considers a phase space distribution $\rho \in {\cal F}(\mathbb{R}\times M)$ whose total/material time derivative
$$ 0~\approx~\frac{d\rho}{dt}~=~\sum_{\mu =0}^{2n}\frac{\partial
\rho}{\partial z^{\mu}} \dot{z}^{\mu} ~\approx~\frac{\partial \rho}{\partial
t}+\{\rho,H\}_{PB}\tag{17} $$
vanishes. Eq. (17) is called Liouville's equation. It expresses conservation of probability. Next we form the current
$$J^{\mu}~=~\rho V^{\mu},\tag{18} $$
which is divergence-free
$$\sum_{\mu=0}^{2n}\frac{\partial J^{\mu}}{\partial z^{\mu}}
~\approx~0\tag{19}$$
because of eqs. (12-18). Eq. (19) can be viewed as a continuity equation. Note that the index $\mu$ runs over dynamically active phase space variables and time. Eq. (19) is quite different from a Noether conservation law
$$\sum_{\mu=0}^{D-1}\frac{d j^{\mu}}{d x^{\mu}} ~\approx~0,\tag{20}$$
where the index $\mu$ runs over dynamically passive spacetime variables. The best bet on a connection with Noether's theorem seems to be eq. (17) and this Phys.SE post.
References:

*

*C. Crnkovic & E. Witten, Covariant description of canonical formalism in geometrical theories. Published in Three hundred years of gravitation (Eds. S. W. Hawking and W. Israel), (1987) 676.

--
$^1$ We use for simplicity here the convention that super-derivatives and the transformation ${\bf s}$ are  left derivation, i.e.
$$ {\bf s}(fg)~=~{\bf s}(f)~g + (-1)^{|f|}f ~{\bf s}(g). \tag{21} $$
A: Regarding your first question, if you pass to complex variables, Noether theorem implies your conservation law. The point is that you are dealing with two independent real solutions while the differential equation refers to only one solution. The simplest way to introduce two independent solutions is viewing them as the real and complex part of a complex solution.
In practice consider the Lagrangian:
$$L = \dot{\overline{z}}\dot{z} + \overline{z}z\tag{1}$$ 
where $z = x_1+ix_2$. The equation of motion is
$$\ddot{z}+ z=0\:,$$
that is
$$\ddot{x}_j+ x_j=0\quad j=1,2\:.$$
The Lagrangian (1) is invariant under the one-parameter group 
$$z \to e^{ia}z\quad a \in \mathbb R\:.$$
Noether theorem gives rise to the complex conserved quantity
$$I = i\dot{\overline{z}}z -  i\dot{z}\overline{z}\:.$$
You easily see that
$$\frac{1}{2}I = x_1\dot{x}_2 - x_2\dot{x}_1\:.$$
