Your model fit and ANOVA table are

    lm = LinearModelFit[data3, {species, sl}, {species, sl}, NominalVariables -> species];
    lm["ANOVATable"]

[![ANOVA Table][1]][1]

The P-value of 9.08206 x 10^-7 looks "suspicious" but that's because that P-value isn't what you think it is.  *Mathematica* is using what is called (in SAS lingo) "Type I Sums of Squares" or "Sequential sum of squares".  This is the sum of squares associated with species unadjusted for the covariate (`sl` = standard length). (And more generally sequential sums of squares are "adjusted for terms preceding it and unadjusted for terms following it".)

You can get an appropriate P-value for testing whether the intercepts for both species are identical by simply changing `{species, sl}, {species, sl}` to ` {sl, species}, {species, sl}`.

    lm = LinearModelFit[data3, {sl, species}, {species, sl}, NominalVariables -> species]
    lm["ANOVATable"]

[![ANOVA table with correct P-values for species equivalence][2]][2]


Does that sound crazy?  Yes, it does.  I love *Mathematica* but it has a long ways to go before I would use it for analysis of variance or analysis of covariance.

**Addition**

If you need to do a lot of multiple comparisons and want to use various adjustments (such as Tukey, Bonferroni, etc.), you should use R or SAS rather than *Mathematica*.  But if you're just interested in all possible pairwise comparisons from a simple ANCOVA, here's one (convoluted) way to do it.

Create some data.

    groups = {"control", "treatment 1", "treatment 2", "treatment 3"};

    data = BlockRandom[SeedRandom[123];
       Block[{vals, times, rand}, vals = RandomChoice[groups, 100];
        times = RandomInteger[10, 100];
        rand = RandomReal[1, 100];
        Transpose[{vals, times, (vals /. 
             Thread[Rule[groups, {.16, .34, .57, 1.1}]]) - .05 times + rand}]]];

Fit the model and get estimates of parameters.

    lm = LinearModelFit[data, {treatment, time}, {treatment, time}, 
      NominalVariables -> treatment, IncludeConstantBasis -> False]
    estimates = lm["BestFitParameters"];
    cov = lm["CovarianceMatrix"];

Now find all possible differences in the treatments (i.e., the intercepts of a simple linear model).

    (* Number of treatments *)
    t = Length[groups]; 

    (* Degrees of freedom for error *)
    df = lm["ANOVATableDegreesOfFreedom"][[Length[ lm["ANOVATableDegreesOfFreedom"]] - 1]]
    (* 95 *)

    (* Indices of all possible pairwise comparisons *)
    comparisons = Select[Sort[#] & /@ Tuples[Range[t], {2}], #[[1]] < #[[2]] &] // DeleteDuplicates
    (* {{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}} *)
    
    (* Differences and standard errors of differences *)
    diffs = estimates[[#[[1]]]] - estimates[[#[[2]]]] & /@ comparisons
    (* {-0.070739, -0.425142, -0.814842, -0.354403, -0.744103, -0.389699} *)

    seDiffs = Sqrt[cov[[#[[1]], #[[1]]]] + cov[[#[[2]], #[[2]]]] - 
         2 cov[[#[[1]], #[[2]]]]] & /@ comparisons
    (* {0.0883031, 0.089601, 0.0895971, 0.0828474, 0.0879348, 0.0891328} *)

Now put together the pieces to determine the Bonferroni adjusted P-values and confidence limits for the pairwise differences:

    tValues = (estimates[[#[[1]]]] - estimates[[#[[2]]]])/
        Sqrt[cov[[#[[1]], #[[1]]]] + cov[[#[[2]], #[[2]]]] - 
          2 cov[[#[[1]], #[[2]]]]] & /@ comparisons
    (* {-0.801093, -4.74484, -9.09451, -4.27779, -8.46198, -4.37212} *)

    Pvalues = 2*CDF[StudentTDistribution[df], -Abs[#]] & /@ tValues
    (* {0.425077, 7.34947*10^-6, 1.42412*10^-14, 0.0000449881,  3.17124*10^-13, 0.0000314736} *)

    α = 0.05/Length[comparisons]
    (* 0.00833333 *)

    confLimits = Transpose[{diffs - InverseCDF[StudentTDistribution[95], 1 - α/2] seDiffs,
       diffs + InverseCDF[StudentTDistribution[95], 1 - α/2] seDiffs}]
    (* {{-0.308683, 0.167205}, {-0.666584, -0.183701}, {-1.05627, 0.573411}, 
        {-0.577646, -0.131161}, {-0.981054, -0.507151}, {-0.629879, -0.14952}} *)



  [1]: https://i.sstatic.net/WTqEN.png
  [2]: https://i.sstatic.net/Dfyba.png