It does follow from calculus. Here's the standard way this is treated (I'm not going to be explicit about mathematical details such as smoothness assumptions here).
Definition of $\delta q$.
Given a parametrized path $q:t\mapsto q(t)$, we consider a deformation of the path which we call $\hat q:(t, \epsilon)\mapsto \hat q(t,\epsilon)$ satisfying $\hat q(t,0) = q(t)$. The parameter $\epsilon$ is the deformation parameter. We can now define the variation $\delta q$ of the path $q$ as follows:
\begin{align}
\delta q(t) = \frac{\partial\hat q}{\partial \epsilon}(t,0) \tag{$\star$}
\end{align}
To motivate this definition, notice that we can Taylor expand $\hat q$ in the $\epsilon$ argument about $\epsilon=0$ as follows:
\begin{align}
\hat q(t,\epsilon) = \hat q(t,0) + \epsilon \frac{\partial\hat q}{\partial\epsilon}(t,0) + O(\epsilon^2)
\end{align}
which, in light of the definition of $\delta q$ above can be rewritten as
\begin{align}
\hat q(t,\epsilon) = q(t) + \epsilon\delta q(t) + O(\epsilon^2)
\end{align}
so that we recognize $\delta q(t)$ as the first order Taylor coefficient of the deformation $\hat q$ when we expand in the deformation parameter. Note that some authors in physics will instead define $\delta q$ with an extra factor of $\epsilon$ on the right hand side of $(\star)$, but this is just a matter of convention.
The commutativity property.
Now that we have defined $\delta q$, we address the commutativity of $\delta$ and $t$-derivatives. Well, now that everything is very explicit, this is pretty straightforward. First, we need note that $\dot q$ is a different curve than $q$, so we need to define its variation $\delta\dot q$. The standard way to do this is to induce this variation using the same deformation $\hat q$. Namely, we define
\begin{align}
\delta\dot q(t) = \frac{\partial^2\hat q}{\partial \epsilon\partial t}(t,0) \tag{$\star\star$}
\end{align}
then we can compute
\begin{align}
\frac{d}{dt}\delta q(t) = \frac{d}{d t}\left(\frac{\partial\hat q}{\partial \epsilon}(t,0)\right) = \frac{\partial^2\hat q}{\partial t\partial \epsilon}(t,0) = \frac{\partial^2\hat q}{\partial \epsilon\partial t}(t,0) = \delta\dot q(t)
\end{align}
which is the desired result.
Naturalness questions.
In some sense, the definitions $(\star)$ and $(\star\star)$ are arbitrary, but only to the extent that any definition is always arbitrary because we have to choose it. They are, however, standard and pretty physical if you ask me.
To get intuition for $(\star)$, consider $\hat q(t,\epsilon)$, and imagine fixing some $t_*$. Then as we vary $\epsilon$, we obtain a curve $\epsilon\mapsto \hat q(t_*, \epsilon)$. The variation $\delta q(t_*)$ is the derivative of this curve with respect to $\epsilon$ evaluated at $\epsilon = 0$, in other words, it is its tangent vector at $\epsilon = 0$ (think velocity). This tangent vector simply tells us the "direction" in which the original curve $q$ is changing at point $t_*$ as we apply the deformation to it. See the following diagram (which I hope is more clear than what I just said)
Here's another way of seeing that the definition $(\star)$ is natural which also shows why $(\star\star)$ is natural. In classical mechanics, we often consider a system described by an action which is the integral of a local lagrangian;
\begin{align}
S[q] = \int dt\,L(q(t), \dot q(t), t).
\end{align}
Now, suppose that we want to determine what happens to $S[q]$ when we deform the path $q$. Using the notation $\hat q$ from above for the deformation, this amounts to evaluating $S[\hat q(\cdot,\epsilon)]$. Let's compute this quantity to first order in epsilon. We find that
\begin{align}
S[\hat q(\cdot, \epsilon)]
&= \int dt\, L\left(\hat q(t,\epsilon), \frac{\partial\hat q}{\partial t}(t,\epsilon), t\right) \\
&= S[q] +\epsilon \int dt\left[\frac{\partial L}{\partial q}(q(t), \dot q(t), t)\delta q(t) + \frac{\partial L}{\partial \dot q}(q(t), \dot q(t), t)\delta \dot q(t)\right] + O(\epsilon^2)
\end{align}
I have skipped some steps here, but the point is that the quantities $\delta q$ and $\delta\dot q$ that we defined in $(\star)$ and $(\star\star)$ naturally arise in the context of taking the variation of a functional of the path $q$. In particular, the variation of $\dot q$ induced by the variation of $q$ as defined in $(\star\star)$ is the object that naturally arises, not some other independent variation.
However, see Qmechanic's answer below which points out that in other contexts, like when using D'Alembert's principle, the variations $q$ and $\dot q$ may not have precisely the same meaning as they do in the contexts described above, and in these contexts the commutativity rule need not hold.
