Context. A variety of formation models for dwarf spheroidal (dSph) galaxies have been proposed in the literature, but generally they have not been quantitatively compared with observations. Aims. We search for chemodynamical patterns in our observational data set and compare the results with mock galaxies consisting of pure random motions, and simulated dwarfs formed via the dissolving star cluster and tidal stirring models. Methods. We made use of a new spectroscopic data set for the Milky Way dSph Leo I, combining 288 stars observed with Magellan/IMACS and existing Keck/DEIMOS data, to provide velocity and metallicity measurements for 953 Leo I member stars. We used a specially developed algorithm called BEACON to detect chemo-kinematical patterns in the observed and simulated data. Results. After analysing the Leo I data, we report the detection of 14 candidate streams of stars that may have originated in disrupted star clusters. The angular momentum vectors of these streams are randomly oriented, consistent with the lack of rotation in Leo I. These results are consistent with the predictions of the dissolving cluster model. In contrast, we find fewer candidate stream signals in mock data sets that lack coherent motions ∼99% of the time. The chemodynamical analysis of the tidal stirring simulation produces streams that share a common orientation of their angular momenta, which is inconsistent with the Leo I data. Conclusions. Even though it is very difficult to distinguish which of the detected streams are real and which are only noise, we can be certain that there are more streams detected in the observational data of Leo I than expected in pure random data.