Quasisymmetry (QS) is a hidden symmetry of the magnetic field strength, B, that confines charged particles effectively in a three-dimensional toroidal plasma equilibrium. Here, we show that QS has a deep connection to the underlying symmetry that makes solitons possible. Our approach uncovers a hidden lower dimensionality of B on a magnetic flux surface, which could make stellarator optimization schemes significantly more efficient. Recent numerical breakthroughs have yielded configurations with excellent volumetric QS and surprisingly low magnetic shear. Our approach elucidates why the magnetic shear is low in these configurations. Furthermore, we deduce an upper bound on the maximum toroidal volume that can be quasisymmetric and verify it for the Landreman-Paul precise quasiaxisymmetric (QA) stellarator configuration. In the neighborhood of the outermost surface, we show that the B approaches the form of the 1-soliton reflectionless potential. We present three independent approaches to demonstrate that quasisymmetric B is described by well-known integrable systems such as the Korteweg-de Vries (KdV) equation. The first approach is weakly nonlinear multiscale perturbation theory, which highlights the crucial role that magnetic shear plays in QS. We show that the overdetermined problem of finding quasisymmetric vacuum fields admits solutions for which the rotational transform is not free but highly constrained. We obtain the KdV equation (and, more specifically, Gardner's equation for certain choices of parameters). Our second approach is non-perturbative and based on ensuring single-valuedness of B, which directly leads to its Painlevé property shared by the KdV equation. Our third approach uses machine learning, trained on a large dataset of numerically optimized quasisymmetric stellarators. We robustly recover the KdV (and Gardner's) equation from the data.