FFT starting at zero
The FFT output will be bins evenly spaced from DC (bin 0) to 1 bin less than your sampling frequency. Thus your sampling frequency is at N for bins 0 to N-1, where N is the number of samples in your time domain sequence. 

Within the first Nyquist zone of $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate, the bins from 0 to N/2-1 are the positive frequency components, while the bins from N/2 to N-1 are the negative components. Specifically bin N-1 is -1, bin N-2 is -2, etc so you could view the domain as a circle from 0 to N-1 and repeating again. 

The reason your peaks are a mirror image is because your time domain signal is real, and for real signals the positive and negative frequencies have a conjugate symmetric relationship (same magnitude, opposite phase).

The reason you see a large response in bin 0 is because your time domain signal has a non-zero mean. Bin 0, being the DC term, is the mean of your signal (scaled by a factor of N). You likely don't care about DC, and would therefore benefit by removing the mean of your signal before taking the FFT. 

The FFT is an algorithm that implements the DFT efficiently, which is as the name implies a discrete version of the Fourier Transform. I think it is helpful to see that the Fourier Transform is simply a correlation as we sweep frequency. (Where correlation is complex conjugate multiply and integrate, or equivalently discrete as complex conjugate multiply and accumulate). Thus for the DFT, each bin is a correlation of a single rotating phasor (your traditional sine wave is made of two equal and opposite rotating phasors as shown by Euler's Identity: $\frac{1}{T}(e^{j\omega t}+e^{-j\omega t}$).  So for the DFT, given as: $$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j(k\omega_0) n}$$

So we see that for the first bin (k=0), we just have N times the mean of x. If x has a DC value, we would get a strong correlation for that calculation. Also observe when k =1, which corresponds to the next bin, the exponential is spinnning at its lowest possible rate $\omega_0$ corresponding to one cycle over your time domain sequence. So in that calculation we are computing the correlation of the lowest possible frequency that is spinning in the opposite direction ($e^{+j\omega t}$): such a signal in our x(t), after this computation would have the same result as bin 0: we would get N times the mean after having "derotated" the signal that was present there. Thus every bin goes through the same "averaging filter" which is a Sinc function in frequency (not a great filter), so our output at any given bin will be sensitive to some degree of frequencies everywhere in our spectrum (and as we see with the Sinc function there will be discrete frequencies that are at null positions and therefore will be locations where there is no "leakage" into other bins). 

This should give you enough insight into some important fundamental points to dig further into the DFT/FFT. In addition here are some other related posts I think will help you:

Frequency measurement from acoustic emission sensor data(VDC RMS)

What proportion of a padded FFT should be actual values

Signal Processing/FFT gives very high magnitudes for low frequencies

You will want to read up a little more about windowing too; the links above and other posts (by search) include more information on that. Ultimately fred harris' paper has been a great resource for me on that topic.