As was mentioned in some of the comments, the Lagrangians
$$
\mathcal{L}=(\partial^\mu\psi)^\dagger(\partial_\mu\psi)-m^2\psi^\dagger\psi
$$
and
$$
\mathcal{L}=(D^\mu\psi)^\dagger(D_\mu\psi)-m^2\psi^\dagger\psi+\frac{1}{4}F_{\mu\nu}F^{\mu\nu}
$$
represent distinct theories each with their own properties.
The usual way to motivate the transition from the "ungauged" theory and the "gauged" one is to note that if we want invariance under the transformation $\psi\rightarrow e^{i\alpha}\psi$ for $\alpha=\alpha(x)$ an arbitrary real function, then taking a Lagrangian that's already invariant in the special case where $\alpha$ is a constant and replacing all the derivatives of $\psi$ by covariant derivatives $D_\mu$, would be good enough to construct a Lagrangian which is also invariant under the local transformations.
There is another way of looking at things, however, which may feel a little less ad hoc. Though this viewpoint can be described in terms of this example of $\psi$ fields, it's slightly more natural to start with the example of a vector field.
So, suppose that $V^a$ are the components of some vector field - note that these are only the components. The vector field itself, meaning the abstract object which is invariant under coordinate changes is $V=V^a\boldsymbol{e}_a$ where the $\boldsymbol{e}_a$ form a basis of vectors at each point in space (technically called frame fields). For example, in two dimensions, we could take $\boldsymbol{e}_0=\boldsymbol{\hat r}$ and $\boldsymbol{e}_1=\boldsymbol{\hat \theta}$.
Now the key assumption is that the physics of our system should not depend on the basis vectors we choose to represent our vector fields in - that is, if we changed to Cartesian unit vectors instead of polar unit vectors, the components $V^a$ would certainly need to change, but the object $V=V^a\boldsymbol{e}_a$ should not.
Since any change in the basis vectors $\boldsymbol{e}_a$ will be a (linear) map from a linear space to itself, these can be represented by matrices $U^a_b$ so under a change in basis we would have $\boldsymbol{e}^\prime_a=U^b_a\boldsymbol{e}_b$. If we are truly to be independent of the basis vectors, we will be able to perform such a transformation point by point, these basis change matrices may have arbitrary dependence on the spacetime point, $U^a_b=U^a_b(x)$. In order for $V$ to be independent of these changes, the components must transform by the inverse of $U$, $V^{\prime\,a}=U^{-1\, a}_b V^b$.
Finally now, we want to build our Lagrangian out of $V$ and its derivatives. So long as our manifold has a metric, we can build arbitrarily high derivatives out of the differential $d$ and the Hodge dual $*$. If we compute the differential of $V$ in terms of the components, we would find
$$
dV=(dV^a)\boldsymbol{e}_b+V^a(d\boldsymbol{e}_b).
$$
The differential of the components is simple because these are all $0$-forms (scalars), and so $d V^a=\partial_\nu V^adx^\nu$. For the differential of the basis vectors, we can first note that the result must
a) be a 1-form
b) be some combination of unit vectors again.
These two statements together imply the differential must take the generic form
$$
d\boldsymbol{e}_a=(A_\mu)_a^b\boldsymbol{e}_bdx^\mu
$$
where $A_{\mu\,b}^a$ is some unknown function, suggestively named. Putting this result back into the calculation of $dV$, we find
$$
dV=\partial_\mu V^a\boldsymbol{e}_adx^\mu+V^aA_{\mu\,a}^b\boldsymbol{e}_bdx^\mu.
$$
Collecting the differentials, unit vectors, and components together, this becomes
$$
dV=\boldsymbol{e}_adx^\mu(\delta^a_b\partial_\mu+A_{\mu\,b}^a)V^b=\boldsymbol{e}_adx^\mu(D_\mu)^a_bV^b.
$$
In the last line we have identified the covariant derivative $D$. This differs slightly from the covariant derivative in the question by overall scalings of $A$ (the $iq$) which could have been absorbed into our definition of $A$.
This expression also differs slightly from what's in the question by the additional indices $a$ and $b$ floating around. In the case of the complex scalar field, we are not dealing with a vector, but instead some object $\tilde \psi=\psi z$ where now $z$ is some complex number with $|z|=1$. This now plays the role our $\boldsymbol{e}$'s played before (but has no indices).
Since $z$ must have modulus 1, we can only transform to a new $z$ by $z^\prime=e^{iq\alpha}z$ where $\alpha=\alpha(x)$ in the same way the change of basis matrix $U$ was allowed to vary point to point (and $q$ has been put in for convenience). Since there are no indices on this $z$, our calculation of the differential would yield
$$
d\tilde \psi=dx^\mu zD_\mu\psi=dx^\mu z(\partial_\mu+iqA_\mu)\psi.
$$
As a fun side note, observe that if in the example of a vector we renamed $A$ to $\Gamma$ and called the gauge potential a Christoffel symbol instead, we would immediately reproduce the covariant derivative from general relativity.