Why can't DFT be used when samples are not equally spaced in time?
You can actually use the DFT on non-unformly sampled data. It is still an orthogonal bijective map and has no inherent mathematical issues. 

However, interpreting the meaning of the resulting spectrum is much less intuitive than for evenly spaced data.

The reason is that the basis functions are no longer eigenfunctions of the time translation operator but of an operator that shifts samples, i.e. each sample gets shifted by a possibly different time offset. This is not a very meaningful operation in general and the eigenvectors of it are also not terribly useful.

As a result the bins of the spectrum do not correspond to a single frequency, but rather a single rate of samples. It is possible to associate a single bin with a sinusoid of time-varying momentary frequency though.

This can still be useful, mostly if your sampling intervals are not random but intentionally spaced differently to follow some kind of underlying structure of your data.

So in total, saying that you can't use a DFT on non uniformly sampled data is not accurate. You can, but you have to know exactly what you're doing and be careful with drawing conclusions.

EDIT:

Let me add some remarks about the relationship between DFT and DTFT in this context. The commend of Dilip in the question you refer to has probably the following background. The discrete time Fourier transform can be written as the continuous time Fourier transform of a series Dirac delta distributions:

$$ \mathrm{DTFT}\{x[n]\} = \mathcal{F} \{ \sum_nx[n]\delta(t-t_n) \}$$

where the $t_n$ are the sampling times associated with the sampled values $x[n]$. This relationship naturally generalises to a nonuniform sampling grid if you just assume that the $t_n$ are not given by $nT$ for a sampling period $T$. I'm not sure I would call the result of this a DTFT anymore, but it's mathematically valid to do. because the properties of the continuous time Fourier transform make sure that the result is well behaved and that the basis functions are still orthogonal and represent a single frequency. This only works if your time axis is all real numbers, because only then the complex exponential functions expanded into are guaranteed to be both complete and orthogonal.

For a finite time axis where the DFT would be applicable, this method does not work. You can theoretically define a map that has all the properties of the DFT on a uniform grid when applied to a non-uniform grid by resampling the data on a uniform grid, applying the ordinary DFT, and considering these two steps as one "nonuniform DFT transform". In this case you will in general loose certain mathematical properties like unitary/orthogonality however.