Shade a ListLinePlot of many series with their deviations

Update: Wrapping all in a function:

ClearAll[deviationslLLP]
deviationslLLP[ave_, dev_, opts: OptionsPattern[]] := Module[{fill = 
  Join @@ (Thread[Range[Length@ave] -> 
   List/@(Length[ave] # + Range[Length@ave])]& /@ {1, 2}), 
  apd = Style[#, Opacity[0]] & /@ (ave + dev), 
  amd = Style[#, Opacity[0]] & /@ (ave - dev)},
 ListLinePlot[Join @@ {ave, apd, amd}, Filling -> fill, opts]]

Examples:

dt = RandomVariate[NormalDistribution[1, .05], {10, 50}];
averages = Mean[# dt] & /@ Range[10];
sdevs = StandardDeviation[# dt] & /@ Range[10];

deviationslLLP[averages, sdevs, Frame -> True, ImageSize -> 500, 
  PlotLegends -> RandomWord["Noun", 10]]]

---

Original answer:

You can wrap ave - dev and ave + dev pieces of your data with Style[#, Opacity[0]]& to avoid using PlotStyle:

dt = RandomVariate[NormalDistribution[1, .1], {10, 50}];
averages = Mean[# dt] & /@ {1, 2, 3};
sdevs = StandardDeviation[# dt] & /@ {1, 2, 3};
aplussd = Style[#, Opacity[0]]&/@(averages + sdevs);
aminussd = Style[#, Opacity[0]]&/@(averages - sdevs);
fillings[n_] := Join @@ (Thread[Range[n] -> List/@(n # + Range[n])]& /@ {1, 2});

ListLinePlot[Join@@{averages, aplussd, aminussd}, Filling -> fillings[Length @ averages]]

If you wish yo use PlotStyle you can to specify three colors, say, PlotStyle->{Red, Green,Purple}] to get

Or PlotStyle->{Automatic, Green,Purple}] to get

Use NormalDistribution[1, .05] instead of NormalDistribution[1, .1] and Range[10] instead of {1,2,3} to get

Add the option PlotStyle -> 63 to get

Maybe use a function like makePlot in MapThread:

makePlot = Function[{avg, dev, col},
             ListLinePlot[{avg, avg - dev, avg + dev},
               PlotStyle -> {col, None, None}, Filling -> {1 -> {2}, 1 -> {3}}]];

plots = MapThread[makePlot, {avgs, devs, {Blue, Orange, Red}}];

Show[plots]

I have the sense that this has been answered before, but a quick look couldn't find it, so here's my quick take on it:

ClearAll[varianceListPlot];
Options[varianceListPlot] =
  Options[ListLinePlot];
varianceListPlot[
    domMainLowHigh :
     {
      Repeated[
       {
        {__?NumericQ},
        {__?NumericQ},
        {__?NumericQ},
        {__?NumericQ}
        }
       ]
      },
    ops : OptionsPattern[]
    ] /; Length[DeleteDuplicates[Map[Length, domMainLowHigh, {2}]]] ==
     1 :=
  With[
   {
    plots =
     Join @@
      Map[
       With[{dom = #[[1]]},
         Map[
          Thread[{dom, #}] &, 
          {#[[2]], #[[2]] + #[[3]], #[[2]] - #[[4]]}
          ]
         ] &,
       domMainLowHigh
       ],
    ps = Replace[OptionValue[PlotStyle], Automatic :> ColorData[97]]
    },
   ListLinePlot[
    plots,
    PlotStyle ->
     Table[
      If[Mod[i, 3] == 1, 
       If[ListQ@ps, 
        ps[[Mod[Floor[i/3] + 1, Length@ps, 1]]],
        Replace[ps[Floor[i/3] + 1], _ps :> ps]
        ],
       None
       ],
      {i, Length@plots}
      ],
    ops,
    Filling ->
     Flatten@
      Table[
       {
        i -> {i + 1},
        i -> {i + 2}
        },
       {i, 1, Length@plots, 3}
       ]
    ]
   ];
varianceListPlot[
   dataDevs : 
    Repeated[{{_?NumericQ, _?
        NumericQ, {Repeated[_?NumericQ, 2]}}, ___}],
   ops : OptionsPattern[]
   ] :=
  varianceListPlot[
   {#[[All, 1]], #[[All, 2]], #[[All, 3, 1]], #[[All, 3, 2]]} & /@ 
    dataDevs,
   ops
   ];
varianceListPlot[
   dd :
    {
     {__?NumericQ},
     {__?NumericQ},
     {__?NumericQ},
     {__?NumericQ}
     },
   ops : OptionsPattern[]
   ] :=
  varianceListPlot[{dd}, ops];
varianceListPlot[
  dataDevs : {{_?NumericQ, _?
      NumericQ, {Repeated[_?NumericQ, 2]}}, ___}, 
  ops : OptionsPattern[]
  ] :=
 varianceListPlot[{dataDevs}, ops]

I provided a two syntaxes; either you pass the domain, main plot, upper variance, and lower variance as a 4-tuple for each data set or you can do it point-wise as {x, y, {above, below}} and it'll unwrap to the prior syntax.

varianceListPlot[
 {
  {
   dom,
   Sin[dom],
   ConstantArray[1, Length@dom],
   ConstantArray[1, Length@dom]
   },
  {
   dom,
   Cos[dom],
   ConstantArray[1, Length@dom],
   ConstantArray[1, Length@dom]
   },
  {
   dom,
   3*Cos[dom/π],
   ConstantArray[1, Length@dom],
   ConstantArray[1, Length@dom]
   }
  }
 ]