Your provided data is very noisy. You can get more information from it if you filter it first.
I will apply a [`LowpassFilter`](http://reference.wolfram.com/language/ref/LowpassFilter.html) and a logarithmic transform on the $y$ values, and scale down the $x$ values. This usually helps the fitting algorithm.

    datat = Transpose[{#[[All, 1]]/1500, 
         Log10[LowpassFilter[#[[All, 2]], .1]]} &@data];
    ListPlot[datat, PlotRange -> All, Joined -> True]

![Mathematica graphics](https://i.sstatic.net/iyUT6.png)

Now, you can perform the multi-peak fitting process from the [linked discussion](http://mathematica.stackexchange.com/questions/26336/how-to-perform-a-multi-peak-fitting).

It's up to you to decide which peaks are signal and which are artifacts. I am providing here the solution with 4 peaks.

    With[{n = 4}, 
     resfunc = 
      peakfunc[A[#], μ[#], σ[#], x] & /@ Range[n] /. 
       model[datat, n][[2]]]
    (* copied from @Silvia's answer with slight modifications *)
    Show@{Plot[Evaluate[resfunc], {x, -5, 10}, 
       PlotStyle -> ({Directive[Dashed, Thick, 
             ColorData["DarkRainbow"][#]]} & /@ 
          Rescale[Range[Length[resfunc]]]), PlotRange -> All, 
       Frame -> True, Axes -> False, ImageSize -> 700], 
      Plot[Evaluate[Total@resfunc], {x, -5, 10}, 
       PlotStyle -> Directive[Thick, Red, Opacity[.5]], PlotRange -> All, 
       Frame -> True, Axes -> False], 
      Graphics[{PointSize[.003], Gray, Line@datat}]}

[![4-peak fitting result][1]][1]

You will of course have to scale the fitted functions back to the original data, but this should be trivial.

  [1]: https://i.sstatic.net/UBq2Z.png