Some users have expressed doubt at the validity of the accepted answer, so let me make it rigorous. To do that, we first need to make sense of the notation $\frac{dy}{dx}$. (Which is not a trivial task.)
We start by changing perspective: instead of thinking of variables $x$ and $y$ as numbers, we think of them as smooth real valued functions on some manifold $M$. So $x\colon M\to \mathbb{R}$ and $y\colon M\to \mathbb{R}$. You might rightly ask: why would we do that and which manifold $M$ are you talking about? The answer to the first question is: because I'd like to talk about the differentials $dy$, $dx$ and say everyday stuff like "$y$ is a constant" or "$y$ is a function of $x$". All of that is impossible inside first order logic + ZFC if we simply interpret $x$ and $y$ as elements of $\mathbb{R}$. Concerning the second question: think of $M$ as the physical state space underlying the problem we are trying to model and $x$ and $y$ as observables. If that sounds too unfamiliar: it's similar to how people in probability theory assume an underlying space of outcomes $\Omega$, in order to talk about random variables. (Which are the things we really care about and historically came before $\Omega$, just like $dx,dy$ historically came before manifolds, but I'm drifting of.)
So, having fixed the background manifold $M$, whenever you hear me say variable, what I mean is a thing of type $M\to \mathbb{R}$.
Definition. Given two variables $x$ and $y$, call $y$ a function of $x$ if $dx$ is not zero almost everywhere and there exists a variable $q\colon M \to \mathbb{R}$ such that
$$
dy = q \cdot dx.
$$
Intuitively, the equation $
dy = q \cdot dx
$ says that the change of $y$ is determined by the change of $x$, i.e. that $y$ depends on $x$.
It's not hard to show that $q$ is uniquely determined by $x$ and $y$, hence we decide to denote it with $\frac{dy}{dx}$ and call it the derivative of $y$ wrt. $x$. It was originally called the differential coefficient, cause that's what it is.
According to this definition, $3^{5x+1}$ is a function of $x$ (assuming $x$ is a true variabel, i.e. $dx\neq 0$), but it's also a function of $5x+1$. On the other hand, $3^{5x+1}$ is not a function of $3$ since $d3=0$ (what we call a constant). In particular $\frac{d3^{5x+1}}{d3}$ is undefined, as user21820 has been pointing out emphatically.
We can now state the chain rule in Leibniz form
Theorem. If $z$ is a function of $y$ and $y$ is a function of $x$, then $z$ is also a function $x$ and their differential coefficients satisfy
$$
\frac{dz}{dx}=\frac{dz}{dy}\cdot \frac{dy}{dx}
$$
The proof of this is trivial.
From this perspective, what the student in my question was trying to do, was to let $z=3^{5x+1}$ and $y=3$. But the theorem does not apply, since $3^{5x+1}$ is not a function $3$. The same point of view can be used in the example discussed here, where a student attempted to differentiate $x^x$ by taking $x$ as inner function. Although that's allowed it just leads to
$$
\frac{dx^x}{dx}=\frac{dx^x}{dx}\cdot \frac{dx}{dx}
$$
which is of not much use.
This is not to say that I don't appreciate the other answers (also users 21820). Taemyr's is just one of the three perspectives that haven been proposed. It might seem like it needs a lot of background to make it rigorous. But consider that mathematicians understood this stuff for at least 200 years without requiring manifolds to formalize it. And consider that the other approaches also require quite some background to make them rigorous (like quantifiers, variable bindings, the idea of dummy/bound variables etc. or derivatives of functions of multiple variables). Each has it advantages and disadvantages and none seems more right than the others, methinks.
