The set of wavefunctions which corresponds to the same state as $\psi$ is just the set of multiples of $\psi$ by a non-zero complex number (respectively, a complex number of absolute value $1$, if you consider only normalized wavefunctions).
The Aharonov-Bohm effect does not change this. What it changes is that the Hilbert space is no longer the space of $L^2$ functions on $M$, but rather the space of $L^2$ sections of a hermitian line bundle, whose connection is the magnetic field.
Concretely, this means that you can only identify a wavefunction with a complex valued function after choosing a gauge on a simply connected patch of $M$. If you do that on two different patches $U $ and $V$, one wavefunction $\psi$ will yield two complex valued functions which differ from a phase on the overlap $U\cap V$ (even though they are still restrictions of the same element of the Hilbert space). This phase is a transition function which is part of defining the quantum system (it is part of the data defining the Hilbert space).
Edit : Wave-functions as section of a fiber bundle
I don't have a specific reference at hand right now, but I would be searching for key-word like "geometric phases", "geometry of the Aharonov-Bohm effect". In geometric quantization, wave-function appear naturally as section of a line bundle (see the discussion here for example), but this is a mathematics thing so it doesn't really address the physical reasoning behind this.
Let me try do give some explanation : in the presence of a magnetic field, the Hamiltonian and the Schrödinger equation are written in terms of the vector potential $\vec A$ :
$$H = \frac{1}{2m}(p-qA)^2$$
To preserve gauge invariance, when we change vector potential $\vec A\to \vec A +\nabla \lambda$, we need to also change the wavefunction by a phase : $\psi \to e^{-iq\lambda }\psi$. Hence : we can only hope for the wave-function to be a complex valued function once a gauge is specified. To see the right framework to describe the wavefunction, we need to look into the geometry of gauge theory.
In the case of the Aharonov-Bohm effect, we can consider the whole of space, with non-vanishing magnetic field inside the flux tube. We can choose a vector potential defined over the whole $\mathbb R^3$ and solve the Schrödinger equation as usual, with wave-functions being just complex-valued functions.
We may also consider the flux tube to be infinitely thin and remove it from our space, which becomes $M = \mathbb R \times (\mathbb R^2 \backslash \{ 0\})$. Here, the magnetic field vanishes uniformly. Again, we can find a gauge defined over the whole $M$ (for example by taking the previous realistic solution and taking the width of the tube to $0$). We cannot choose $A = 0$ however, because the integral of $A$ along a path encircling the flux tube is the magnetic flux $\Phi\neq 0$.
We can also divide $M$ into two simply connected (open) subsets $U$ and $V$ (the left and right side of the tube). On both $U$ and $V$, we can take $A = 0$ (which makes solving the Schrödinger equation easier !). Formally, we do this by starting with the previous vector potential and performing a gauge transformation with $\nabla \lambda_{U,V} = -A$ on $U$ (resp. on $V$). The transformed wavefunctions are $\psi_{U,V} = e^{-iq\lambda_{U,V}}\psi$. On the overlap $U\cap V$, those two wavefunctions satisfy $\psi_V = e^{iq(\lambda_U-\lambda_V)} \psi_U$.
Mathematically, this is the formula relating two different trivialization of a hermitian line bundle : at each point on $M$, the wave-function takes it value on a different copy of $\mathbb C$; the vector potential is the connection which allows us to differentiate the wave function; on a small enough open subsets, you can identify the wavefunctions with complex valued functions, but you have to change gauge and multiply by a phase when switching from one of those small patch to another.