I feel that framing this question in terms of matrices rather than bilinear forms on a vector space obscures what is actually going on and makes it harder to understand what needs to be proved. Here is how I would describe the problem and the partial answer that results from this description:
Let $V$ be a finite dimensional vector space over a field $k$, and let $\mathsf{B}(V)$ denote the vector space over $k$ consisting of bilinear forms on $V$, i.e., an element $\beta\in \mathsf{B}(V)$ is a bilinear mapping $\beta: V\times V\to k$. An element $\alpha\in \mathsf{B}(V)$ is said to be alternating if $\alpha(x,x) = 0$ for all $x\in V$, and an element $\sigma\in \mathsf{B}(V)$ is said to be symmetric if $\sigma(x,y)=\sigma(y,x)$ for all $x,y\in V$. The subset $\mathsf{A}(V)\subset\mathsf{B}(V)$ of alternating forms is a subspace, as is the subset $\mathsf{S}(V)\subset\mathsf{B}(V)$ of symmetric forms. When the characteristic of $k$ is not $2$, there is a $GL(V)$-invariant direct sum decomposition $\mathsf{B}(V) = \mathsf{A}(V)\oplus\mathsf{S}(V)$. When the characteristic of $k$ is $2$, one has, instead, $\mathsf{A}(V)\subset \mathsf{S}(V)\subset\mathsf{B}(V)$, and, apparently, these inclusions have no $GL(V)$-invariant splittings.
An element $\beta\in \mathsf{B}(V)$ is nondegenerate if, for each $x\not=0$ in $V$, there exists a $y\in V$ such that $\beta(x,y)\not=0$. If $\alpha\in \mathsf{A}(V)$ is nondegenerate, then the dimension of $V$ over $k$ must be even. Conversely, if the dimension of $V$ over $k$ is even, then there exists a nondegenerate $\alpha\in\mathsf{A}(V)$, and, moreover, any other nondegenerate $\overline\alpha\in\mathsf{A}(V)$ is of the form $\overline\alpha = m^\ast(\alpha)$ for some $m\in GL(V)$, where, by definition,
$$
\bigl(m^\ast(\alpha)\bigr)(x,y) := \alpha(mx,my)
$$
for any $m:V\to V$. When $\dim_k(V) = 2n$, let $\mathsf{K}(V) = \Lambda^{2n}(V^\ast)$ denote the $1$-dimensional vector space consisting of $2n$-multilinear alternating functions on $V$. There exists a canonical polynomial mapping $\mathrm{Pf}:\mathsf{A}(V)\to \mathsf{K}(V)$ of degree $n$ that satisfies $\alpha^n = n!\ \mathrm{Pf}(\alpha)$ (a property that defines $\mathrm{Pf}$ when $n$ is less than the characteristic of $k$). This Pfaffian vanishes if and only if $\alpha$ is degenerate, and it satisfies $\mathrm{Pf}\bigl(m^\ast(\alpha)\bigr) = \det(m)\ \mathrm{Pf}(\alpha)$.
When $\alpha$ is nondegenerate, let $Sp(\alpha)\subset GL(V)$ denote the subgroup consisting of those $m\in GL(V)$ such that $\alpha(mx,my)=\alpha(x,y)$ for all $x,y\in V$. Define two subspaces ${\frak{s}}(\alpha)\subset \mathrm{End}(V) \simeq V\otimes V^\ast$ and ${\frak{a}}(\alpha)\subset \mathrm{End}(V)$, by saying that $s\in{\frak{s}}(\alpha)$ if $s^\flat(x,y) := \alpha(x,sy)$ is symmetric, while $a\in{\frak{a}}(\alpha)$ if $a^\flat(x,y) := \alpha(x,ay)$ is alternating. Note that ${\frak{s}}(\alpha)$ is a subalgebra of $V\otimes V^*$ under the commutator bracket; in fact, it is the Lie algebra of $Sp(\alpha)$. The subspaces ${\frak{s}}(\alpha)$ and ${\frak{a}}(\alpha)$ are invariant under conjugation by elements of $Sp(\alpha)$.
Now, there is a natural map $S:{\frak{s}}(\alpha)\to {\frak{a}}(\alpha)$, given by $S(s) = s^2$. In other words, if $\alpha(x,sy)$ is symmetric, then $\alpha(x,s^2x) = \alpha(sx,sx) = 0$, so $\alpha(x,s^2x)$ is alternating.
Here, then, is the question: What is the image of $S$? (The OP is actually asking whether the image of $S$ contains the invertible elements of ${\frak{a}}(\alpha)$.)
Note that the dimension of ${\frak{s}}(\alpha)$ is $2n^2{+}n$, while the dimension of ${\frak{a}}(\alpha)$ is $2n^2{-}n$, so it's conceivable that $S$ is actually surjective.
Also, the map $S$ is $Sp(\alpha)$-equivariant, so it's a question that can be studied by looking at the orbits of this group acting on ${\frak{a}}(\alpha)$.
Remark: It took me a while to recognize that this is what is going on because the question, as asked, sneaks in an extraneous quadratic form that breaks the symplectic symmetry. A 'reference' alternating form on $k^{2n}$ has been specified by the formula $\alpha_0(x,y) = x^TJy$ for $x,y\in k^{2n}$. Note that the matrix $J$ satisfies $J^2 = -I$, an identity that has no meaning for an alternating form. The only way one can interpret an alternating form as a linear transformation (so that squaring makes sense) is to have some other way of identifying $V$ with $V^\ast$. Of course, this is supplied by the linear map $x\mapsto x^T$ in the formula. In other words, a (symmetric) bilinear form $\beta(x,y) = x^Ty$ has been introduced into the picture, and it breaks the symplectic symmetry. Anyway, writing $\alpha(x,y) = x^TAy$ and asking whether one can write $A = MJM$ for $M$ symmetric can be re-interpreted as follows: Note that $M=Js$ where $s\in{\frak{s}}(\alpha_0)$ and that $\alpha(x,y) = x^TAy = x^TJJ^{-1}Ay = \alpha_0(x,J^{-1}Ay) = \alpha_0(x,ay)$ where $a = J^{-1}A$ lies in ${\frak{a}}(\alpha_0)$. Putting this together says that
$$
Ja = A = MJM = (Js)J(Js) = -Js^2,
$$
so showing that the equation $A = MJM$ can be solved is equivalent to showing that the equation $a = -s^2$ for a given $a\in{\frak{a}}(\alpha_0)$ can be solved for some $s\in{\frak{s}}(\alpha_0)$. (It's off by a minus sign, but that's OK because the goal is to characterize the image of $S$, so characterizing its negative is just as good.)
Anyway, back to the question: One approach is to look at the orbits of $Sp(\alpha)$ acting on ${\frak{s}}(\alpha)$ and see what their squares look like. This may be easier because the adjoint orbits of $Sp(\alpha)$ on its Lie algebra have been much studied.
As an example, suppose that $k$ is algebraically closed and (for my comfort) that it has characteristic zero. Say that a pair of nondegenerate alternating $2$-forms $(\alpha_0,\alpha)$ on $k^{2n}$ is generic if the $n$ roots of the equation $\textrm{Pf}(\alpha - \lambda\ \alpha_0) = 0$ are all distinct. Then one can prove (see below) that a basis of $1$-forms on $k^{2n}$ exists so that
$$
\alpha_0 =\theta^1\wedge\theta^2+\theta^3\wedge\theta^4+\cdots
+\theta^{2n-1}\wedge\theta^{2n}
$$
while
$$
\alpha =\lambda_1\ \theta^1\wedge\theta^2+\lambda_2\ \theta^3\wedge\theta^4+\cdots
+\lambda_n\ \theta^{2n-1}\wedge\theta^{2n}.
$$
Thus, the problem uncouples into $n$ separate problems that are each trivially solvable, so, the problem is solvable for the generic pair in this case.
As another example, in the case $n=2$, for an arbitrary field (even one of characteristic $2$), one can, by hand, classify the pairs $(\alpha_0,\alpha)$ with $\alpha_0$ nondegenerate and show that $S:{\frak{s}}(\alpha_0)\to{\frak{a}}(\alpha_0)$ is surjective. (I'll put in the details if someone asks.)
To prove surjectivity for all $n$, one may need to understand the orbits of $Sp(\alpha)$ acting on ${\frak{a}}(\alpha)$. I think that this is a classical problem (I'm not an algebraist, so I'm not completely sure), so maybe it's time to look at the literature. The classification of the possible $Sp(\alpha)$-orbit types in ${\frak{a}}(\alpha)$ gets more complicated as $n$ increases, so maybe some other approach needs to be tried. One would expect the orbits of $Sp(\alpha)$ in ${\frak{a}}(\alpha)$ to be somewhat simpler than the orbits of $Sp(\alpha)$ in ${\frak{s}}(\alpha)$, just because the dimension is lower. However, I note that the rings of $Sp(\alpha)$-invariant polynomials on each of the vector spaces ${\frak{s}}(\alpha)$ and ${\frak{a}}(\alpha)$ are each free polynomial rings on $n$ generators, so it may be that the complexity of the orbit structures are (at least roughly) comparable in the two cases.
The uncoupling step: In the general case, for $a\in{\frak{a}}(\alpha_0)$, one has $\textrm{Pf}(a^\flat - \lambda\ \alpha_0) = p_a(\lambda)\ \textrm{Pf}(\alpha_0)$, where $\det(a - \lambda I) = p_a(\lambda)^2$. Letting $f_1(\lambda),\ldots,f_k(\lambda)$ denote the distinct irreducible factors of $p_a(\lambda)$, one has
$$
p_a(\lambda) = f_1(\lambda)^{d_1}\cdots f_k(\lambda)^{d_k}.
$$
There is a direct sum decomposition $V = V_1\oplus\cdots\oplus V_k$ into the corresponding generalized eigenspaces of $a$, and one sees without difficulty (using the identity $\alpha_0(ax,y) = \alpha_0(x,ay)$) that $\alpha_0(V_i,V_j) = 0$ for $i\not=j$. Moreover, any solution $s\in{\frak{s}}(\alpha_0)$ of $a = \pm s^2$ must commute with $a$ and therefore preserve its generalized eigenspaces. Thus, generalizing the 'uncoupled' situation given in the first example, one sees that the problem reduces to the case in which $p_a(\lambda)$ is a power of a single irreducible polynomial.
Unfortunately, it turns out that, even in this uncoupled case, the minimal polynomial of $a$ can fail to be irreducible (i.e., $a$ need not be semi-simple), and I do not know a simple way to handle all of these cases. When $n=2$, this can be handled 'by hand', but even for $n=3$, it seems to be a little tricky (although, I think that I have correctly handled the cases there and shown surjectivity in that case as well).
The uncoupled semi-simple case: Here is how one can complete the proof of solvability in the uncoupled, semi-simple case, i.e., when the minimal polynomial of $a$ is irreducible. Let $p(\lambda)=0$ be the (irreducible) minimal polynomial of $a$, say of degree $m$. We can assume that $p(0)\not=0$, since, otherwise, $a=0$, and there is nothing to prove (i.e., one can just take $s=0$, and the problem is solved). Let $K\subset End(V)$ denote the field generated by $a$, so that $[K:k]=m$.
Now, for any nonzero $x\in V$, one has $\alpha_0(K{\cdot}x, K{\cdot}x) = 0$, since $\alpha_0(a^ix,a^jx)=0$ for all $i$ and $j$. If one chooses $x,y\in V$ such that $\alpha_0(x,y)\not=0$, then it is easy to see (using $p(0)\not=0$) that $\alpha_0$ is nondegenerate on the $a$-invariant subspace $W = K{\cdot}x\oplus K{\cdot}y$. One then can write $V = W\oplus W^\perp$ (where the $\perp$ is taken with respect to $\alpha_0$) and see that the problem uncouples into separate problems on $W$ and $W^\perp$. By induction, it then suffices to consider the case $W = V$ (and, hence, $n=m$).
In this case, the general element of $V = W = K{\cdot}x\oplus K{\cdot}y$ can be written uniquely in the form $z = f_1(a)x + f_2(a)y$ where $f_1$ and $f_2$ are polynomials (with coefficients in $k$) of degree at most $m{-}1$. Let $q(a)$ and $r(a)$ be any polynomials in $a$ (i.e., elements of $K$) and define a map $s:V\to V$ by
$$
s\bigl(f_1(a)x + f_2(a)y\bigr) = f_1(a)q(a)y + f_2(a)r(a)x.
$$
One checks that $s\in{\frak{s}}(\alpha_0)$ and notes that, by construction, one has the identity $s^2 = q(a)r(a)$. By setting $q(a) = a$ and $r(a) = \pm1$, one sees that one may arrange $s^2 = \pm a$. This $s$ solves the problem.
