There are already some good answers here, but I just wanted to add an approach to thinking about natural variables using the notion of Availability- which I think makes the whole thing more intuitive.
Derivation of the Availability as extractible work
We start from the entropy of a system, and the principle that entropy must always increase:
$$ dS \ge 0 $$
For an isolated system this expression is useful enough by itself, and we can just rearrange the first law to get:
$$ dS = \frac{dU + pdV - \mu dN}{T} $$
The equation of state can then be subbed into this equation, which can then be solved when the various constraints are applied.
When there is a system in contact with a reservoir, which is often how systems are thought of in thermodynamics, it is much more convenient to consider the change of entropy of the system and the reservoir simultaneously. To do this we can define the Availability of a system, using the second law, that entropy must always increase as stated above, as:
$$ dS_{tot} = dS_{sys} + dS_{reservoir} \ge 0 $$
$$ dS_{tot} = dS_{sys} + \frac{dU_{R} + p_{R}dV_{R} - \mu_{R} dN_{R}}{T_{R}} \ge 0 $$
$$ dS_{tot} = \frac{dS_{sys}T_{R} + dU_{R} + p_{R}dV_{R} - \mu_{R} dN_{R}}{T_{R}} \ge 0 $$
Using the fact that the system and reservoir are isolated from any outside systems, we have $ dV_{R} = - dV_{sys}$, $ dN_{R} = - dN_{sys}$ and $ dU_{R} = - dU_{sys}$. Subbing these in then simplifies the above expression to:
$$ dS_{tot} = \frac{dS_{sys}T_{R} - dU_{sys} - p_{R}dV_{sys} - \mu_{R} dN_{sys}}{T_{R}} \ge 0 $$
Now, dropping the subscript for the 'system' variables, the availability can be written as:
$$ dA = - T_{R}dS_{tot} = dU + p_{R}dV + \mu_{R} dN - T_{R}{dS} $$
$$ dA = (T -T_{R})dS - (p -p_{R})dV + (\mu -\mu_{R})dN $$
Where all the changes on the right hand side are with respect to the system variables.
The second law still applies to the availability, but in this new form it is written as:
$$ dA \le 0 $$
Now all this might seem superfluous to an explanation of the natural variables so far, but if we consider what the availability actually means then we can distill the intuition behind this reconstruction of the entropy.
Consider the amount of work which can be extracted from a system and reservoir ensemble. The amount of extractable work is clearly just the sum of internal energy in the system and reservoir:
$$ dW_{extractable} = -dU - dU_{R} $$
i.e. an amount of work done by the ensemble must be equal to the change in internal energy in the system and the reservoir- just conservation of energy. Now if we sub in the availability expression derived above we have:
$$ dW_{extractable} = p_{R}dV + \mu_{R} dN - T_{R}{dS} - dA - dU_{R} $$
Subbing in $ dV_{R} = - dV$, $ dN_{R} = - dN$, and the definition of dU from the first law, and rearranging gives:
$$ dW_{extractable} = -dA - T_{R}d(S + S_{R}) $$
In the reversible case the second term on the right handside is equal to zero, and so we have:
$$ - dA = dW_{extractable, max} $$
The availability gives the maximum amount of extractable work from a system.
Connection to Natural Variables
The amount of extractable work is of primary importance in many situations, though it often not explicitly stated in terms of the Availability. This is because historically, thinking about the change in free energy of an experimental system is much more convenient, e.g. Gibbs Free Energy is often used in chemistry, and the Hemholtz Free Energy in Physics. When these terms are introduced to students they are often just defined without reference to their derivations- at least in my experience of school and university Physics and Chemistry. But what is the difference between all these different uses of the term 'Free Energy'?
Essentially, these terms, as well as the Enthalpy and Grand Potential all express the same thing as the Availability, i.e. the amount of extractible work, but with respect to different constraints on the system. And these constraints are what define the natural variables of a state variable. Or more precisely, the natural variables will define the magnitude of a state variable, while the non-natural variables will define the magnitude of a change in the state variable.
E.g. for the Hemholtz free energy we can consider the amount of extractable work at constant T, V and N:
$$ dA = (dU + p_{R}dV + \mu_{R} dN - T_{R}{dS})_{const T, V, N} $$
$$ dA = d(U -TS) $$
$$ A_{const T, V, N} = F = U - TS $$
Where the constant terms are dropped between the first and second line. Thus we have $ F(T, V, N) = U - TS $ and $ dF = -SdT - pdV + \mu dN $
Similarly for the Gibbs free energy, which is the extractable work at constant T, p, N:
$$ dA = (dU + p_{R}dV + \mu_{R} dN - T_{R}{dS})_{const T, p, N} $$
$$ dA = d(U -TS +pV) $$
$$ A_{const T, p, N} = G = U - TS + pV $$
Which gives us $ dG(T, p, N) = \mu dN - SdT + Vdp $.
Finally, the internal energy is the amount of extractable work at constant V, N and S:
$$ dA = (dU + p_{R}dV + \mu_{R} dN - T_{R}{dS})_{const V, S, N} $$
$$ dA = d(U) $$
$$ A_{const V, S, N} = U $$
Which gives us $ dU = TdS -pdV +\mu dN $, which comes from the first law.