Harmonics of square wave
Here is one way to think about this:


Start with the rectangular pulse of 10ms length. DFT is the sinc() function with a "period" of 100Hz
Now repeat periodically with 100ms intervals. Mathematically that's a convolution with a 100ms pulse train.
The DFT of the pulse train is also a pulse train in the frequency domain with 10 Hz spacing
Convolution with in the time domain is multiplication in the frequency domain. Hence we need to multiply the 100 Hz sinc() with a 10 Hz pulse train. That's the basic shape you have: a line spectrum with 10 Hz line spacing with the overall contour of a 100Hz sinc()
Sample at 1kHz. The sinc() function decays fairly slowly with frequency: around the Nyquist frequency you are only at -25 dB or so. Hence you will get non-trivial aliasing. However, since your sample rate is an integer multiple of the line spacing, the aliasing products will fall on top of the existing lines and will be hard to see. It'll just be a small deviation from the sinc() contour.


The "fuzz" between the lines appears to be noise. Your time domain wave form looks noisy and not particularly clean. It could also be caused by the way you calculate the spectrum (windowing, periodicity misalignment, FFT size, etc.)