Dedicated to Fa Yuen Wu on the occasion of his 80th birthday.Abstract. Focusing on examples associated with holonomic functions, we try to bring new ideas on how to look at phase transitions, for which the critical manifolds are not points but curves depending on a spectral variable, or, even, fill higher dimensional submanifolds. Lattice statistical mechanics, often provides a natural (holonomic) framework to perform singularity analysis with several complex variables that would, in the most general mathematical framework, be too complex, or simply could not be defined. In a learn-by-example approach, considering several Picard-Fuchs systems of two-variables "above" Calabi-Yau ODEs, associated with double hypergeometric series, we show that D-finite (holonomic) functions are actually a good framework for actually finding properly the singular manifolds. The singular manifolds are found to be genus-zero curves. We, then, analyse the singular algebraic varieties of quite important holonomic functions of lattice statistical mechanics, the n-fold integrals χ (n) , corresponding to the n-particle decomposition of the magnetic susceptibility of the anisotropic square Ising model. In this anisotropic case, we revisit a set of so-called "Nickelian singularities" that turns out to be a two-parameter family of elliptic curves. We then find a first set of non-Nickelian singularities for χ (3) and χ (4) , that also turns out to be rational or ellipic curves. We underline the fact that these singular curves depend on the anisotropy of the Ising model, or, equivalently, that they depend on the spectral parameter of model. This has important consequences on the physical nature of the anisotropic χ (n) 's which appear to be highly composite objects. We address, from a birational viewpoint, the emergence of families of elliptic curves, and of Calabi-Yau manifolds on such problems. We also address the question of the singularities of non-holonomic functions with a discussion on the accumulation of these singular curves for the non-holonomic anisotropic full susceptibility χ.