diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8d40c73..d838cdf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -22,6 +22,15 @@ first entry is written.
   in-place update over `pl.x` and `pl.v` under `--track-allocation`. It measures rather than
   asserts, so `test/` was the wrong home: a profiling harness that runs as part of the suite costs
   a minute and checks nothing.
+- **`ParticleList` carries a species table and a per-particle species index.** A new `Species`
+  type holds a name, a charge and a mass, and rejects a mass that is not positive (`NaN`
+  included). `charge(pl, i)`, `mass(pl, i)` and `species(pl, i)`
+  read them back for each particle, and `chargeweights(pl)` returns the charge-weighted weights
+  `q_i · w_i` as a lazy `1 × np` matrix, the weight a charge deposit reads. A list built without
+  a species holds one species with charge and mass one, so `q · w` is `w` and every existing
+  result is unchanged. `h5save` writes the species table and the index as a group beside
+  `particles`, and a file without that group loads with the default species. The species index
+  costs 8 bytes per particle at construction.
 
 ### Changed
 
@@ -49,7 +58,8 @@ first entry is written.
   `eachslice` over each variable's view, yielding a zero-allocation `AbstractVector` with
   the same element type and values. Measured at `np = 10⁶` particles, `nd = 6`
   components, three variables: allocation reduced from 328.0 MB to 184.0 MB (43.9%
-  reduction).
+  reduction). The species index, which the constructor builds as well, adds 8 bytes per
+  particle on top of that; see the species entry under New Features.
 - **`StaticArrays` removed from `[deps]` in Project.toml.** The `using StaticArrays: MVector` import
   in `src/ParticleMethods.jl` was the only use in `src/` and was removed by the `Particle(DT, len)`
   change. It remains in `[compat]` and as a test-only dependency in `[extras]` because the test suite
@@ -76,6 +86,11 @@ first entry is written.
   test dependencies are in `test/Project.toml`; `Random` is added, and the tests use
   `SafeTestsets`. The root `Project.toml` has no
   `[extras]`, `[targets]` or test-only compat entries. The suite keeps its 81 tests.
+- **HDF5 floor raised to 0.16.16** (`[compat]` `"0.16.16, 0.17"`). `h5save` now writes the species
+  names as a variable-length string dataset (`species/name`). HDF5.jl 0.16.0 to 0.16.15 cannot read
+  a variable-length string on Julia 1.10, because its loader calls `Base._memcpy!`, which Julia 1.10
+  does not define; 0.16.16 guards that call. A Julia 1.10 environment with an older HDF5 resolves to
+  0.16.16 or later.
 
 ### Bug Fixes
 
@@ -106,3 +121,7 @@ first entry is written.
   (`AbstractVector`), but this constructor's concrete return type is now correct. This changes
   equality semantics: `Particle(Float64, 7) == Particle(MVector{7}(zeros(7)))` was `true` before
   and is now `false`, because `==` requires the state types to match.
+- **`species` and `speciesindex` are field names of `ParticleList`, so they are reserved.** The
+  constructor asserts that variable and parameter names are disjoint from the field names, so a
+  variable or parameter named `species` or `speciesindex` now fails with an `AssertionError`. Rename
+  such a variable or parameter.
diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md
new file mode 100644
index 0000000..af238da
--- /dev/null
+++ b/KNOWN_ISSUES.md
@@ -0,0 +1,13 @@
+# Known issues
+
+### K1 · `ChargeWeights` and `chargeweights` are exported API with no caller.
+
+- **location:** `src/particle_list.jl:199`
+- **evidence:** No code in this package or in its dependents calls `chargeweights`. The type and
+  the function add about 40 lines of source and about 46 lines of tests. A deposit that weights
+  each particle by its charge can write `charge(pl, i) * pl.list[wrow, i]` inline: at `np = 10⁶`
+  the inline form and `chargeweights(pl)[1, i]` measured 692 µs and 695 µs, in one run. Removing
+  the pair after a release is a breaking change, so the decision is cheapest before the next
+  release: remove it, or keep it until the first deposit that uses a species charge calls it.
+- **kind:** dead code
+- **found:** PR #20
diff --git a/Project.toml b/Project.toml
index c8eb66d..7413bfd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,5 +7,5 @@ version = "0.1.2"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 
 [compat]
-HDF5 = "0.16, 0.17"
+HDF5 = "0.16.16, 0.17"
 julia = "1.10"
diff --git a/src/ParticleMethods.jl b/src/ParticleMethods.jl
index 67b86ca..c810aec 100644
--- a/src/ParticleMethods.jl
+++ b/src/ParticleMethods.jl
@@ -11,6 +11,13 @@ using HDF5: H5DataStore
 
 include("hdf5_utils.jl")
 
+include("species.jl")
+
+export Species
+export charge
+export mass
+export species
+
 include("particle.jl")
 
 export Particle
@@ -19,5 +26,6 @@ include("particle_list.jl")
 
 export ParticleList
 export eachparticle
+export chargeweights
 
 end # module
diff --git a/src/particle_list.jl b/src/particle_list.jl
index 50ef860..303d2bd 100644
--- a/src/particle_list.jl
+++ b/src/particle_list.jl
@@ -8,9 +8,12 @@ _sort_names(::NamedTuple{N, T}) where {N, T} = Tuple(sort([N...]))
 _sort_ntuple(nt::NamedTuple) = NamedTuple{_sort_names(nt)}(nt)
 
 """
-    ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple())
-    ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple())
-    ParticleList(x, v, w; parameters = NamedTuple())
+    ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(),
+                 species = [Species(:default, 1, 1)], speciesindex = nothing)
+    ParticleList(DT, np, nd; variables = NamedTuple(), parameters = NamedTuple(),
+                 species = [Species(:default, 1, 1)], speciesindex = nothing)
+    ParticleList(x, v, w; parameters = NamedTuple(), species = [Species(:default, 1, 1)],
+                 speciesindex = nothing)
     ParticleList(h5::H5DataStore, path = "/")
     ParticleList(fpath::AbstractString, path = "/")
 
@@ -20,7 +23,15 @@ particle, and `pl[i]` is that particle as a [`Particle`](@ref).
 `variables` maps names to row indices or row ranges of `list`, and `parameters` maps names to
 values. Each name is a property of the list: `pl.x` is a view into the rows of a variable `x` for
 all particles, and the stored value for a parameter. The names must not repeat each other or the
-field names `list`, `views`, `params`, `particles`, `variables` and `indices`.
+field names `list`, `views`, `params`, `particles`, `variables`, `indices`, `species` and
+`speciesindex`.
+
+`species` is the table of particle species, each a [`Species`](@ref) with a charge and a mass, and
+`speciesindex` gives the species of each particle as an index into that table. The table is
+converted to the element type of `list`, and the index is copied. A list built without `species`
+holds one species with charge and mass one, so `chargeweights(pl)[1, i]` is the weight `w_i` of
+particle `i`. A list with more than one species needs an explicit `speciesindex` of length `np`,
+with every index in `1:length(species)`.
 
 The second form makes a zero list of element type `DT` with `np` particles of `nd` components.
 The third form stacks positions `x`, velocities `v` and weights `w`, and defines the variables
@@ -35,24 +46,44 @@ struct ParticleList{T, ST <: AbstractMatrix{T}, VT <: NamedTuple, PT <: NamedTup
     particles::PART
     variables::VART
     indices::IND
+    species::Vector{Species{T}}
+    speciesindex::Vector{Int}
 
-    function ParticleList(list, views, params, particles, variables, indices)
+    function ParticleList(list, views, params, particles, variables, indices, species,
+            speciesindex)
         @assert isdisjoint(fieldnames(ParticleList), keys(views))
         @assert isdisjoint(fieldnames(ParticleList), keys(params))
         @assert isdisjoint(keys(views), keys(params))
+        length(speciesindex) == size(list, 2) || throw(ArgumentError(
+            "speciesindex has length $(length(speciesindex)), expected $(size(list, 2))"))
+        all(in(eachindex(species)), speciesindex) || throw(ArgumentError(
+            "speciesindex lies outside 1:$(length(species))"))
         new{eltype(list), typeof(list), typeof(views), typeof(params),
             typeof(particles), typeof(variables), typeof(indices)}(
-            list, views, params, particles, variables, indices)
+            list, views, params, particles, variables, indices, species, speciesindex)
     end
 end
 
-function ParticleList(list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple())
+function ParticleList(
+        list::AbstractMatrix; variables = NamedTuple(), parameters = NamedTuple(),
+        species = [Species(:default, one(eltype(list)), one(eltype(list)))],
+        speciesindex = nothing)
     svariables = _sort_ntuple(variables)
     views = map(idx_range -> view(list, idx_range, :), svariables)
     vars = map(v -> eachslice(v; dims = ndims(v)), views)
     particles = [Particle(p; variables = svariables, parameters = parameters)
                  for p in eachcol(list)]
-    ParticleList(list, views, parameters, particles, vars, svariables)
+
+    if speciesindex === nothing
+        length(species) > 1 && throw(ArgumentError(
+            "a ParticleList with $(length(species)) species needs an explicit speciesindex"))
+        index = ones(Int, size(list, 2))
+    else
+        index = collect(Int, speciesindex)
+    end
+
+    ParticleList(list, views, parameters, particles, vars, svariables,
+        Species{eltype(list)}.(species), index)
 end
 
 function ParticleList(DT::DataType, np::Int, nd::Int; kwargs...)
@@ -84,7 +115,9 @@ function Base.:(==)(pl1::ParticleList{T1, ST1}, pl2::ParticleList{
         && pl1.params == pl2.params
         && pl1.particles == pl2.particles
         && pl1.variables == pl2.variables
-        && pl1.indices == pl2.indices)
+        && pl1.indices == pl2.indices
+        && pl1.species == pl2.species
+        && pl1.speciesindex == pl2.speciesindex)
 end
 
 @inline function Base.hasproperty(::ParticleList{T, ST, VT, PT}, s::Symbol) where {
@@ -135,6 +168,61 @@ eachparticle(pl::ParticleList) = pl.particles
 
 Base.iterate(pl::ParticleList, state...) = iterate(pl.particles, state...)
 
+"""
+    species(pl::ParticleList, i::Integer)
+
+Return the [`Species`](@ref) of particle `i` in `pl`, that is `pl.species[pl.speciesindex[i]]`.
+"""
+@inline species(pl::ParticleList, i::Integer) = pl.species[pl.speciesindex[i]]
+
+"""
+    charge(pl::ParticleList, i::Integer)
+
+Return the charge of particle `i` in `pl`, the charge of its [`Species`](@ref).
+"""
+@inline charge(pl::ParticleList, i::Integer) = charge(species(pl, i))
+
+"""
+    mass(pl::ParticleList, i::Integer)
+
+Return the mass of particle `i` in `pl`, the mass of its [`Species`](@ref).
+"""
+@inline mass(pl::ParticleList, i::Integer) = mass(species(pl, i))
+
+"""
+    ChargeWeights{T} <: AbstractMatrix{T}
+
+The charge-weighted weights `q_i · w_i` of the particles in a [`ParticleList`](@ref), as a lazy
+`1 × np` matrix whose `[1, i]` element is the charge of particle `i` times its weight. Returned by
+[`chargeweights`](@ref).
+"""
+struct ChargeWeights{T, PL <: ParticleList} <: AbstractMatrix{T}
+    parent::PL
+    row::Int
+end
+
+Base.size(cw::ChargeWeights) = (1, length(cw.parent))
+
+@inline function Base.getindex(cw::ChargeWeights, i::Int, j::Int)
+    @boundscheck checkbounds(cw, i, j)
+    return charge(cw.parent, j) * cw.parent.list[cw.row, j]
+end
+
+"""
+    chargeweights(pl::ParticleList)
+
+Return the charge-weighted weights of the particles in `pl` as a lazy `1 × np` matrix, whose
+`[1, i]` element is `charge(pl, i) * w_i`, the weight of particle `i`. This is the weight a charge
+deposit reads, so a deposit that weights each particle by its charge changes one argument and no
+signature.
+
+It throws an `ArgumentError` when `pl` has no variable `w`, or when `w` spans more than one row.
+"""
+function chargeweights(pl::ParticleList{T}) where {T}
+    haskey(pl.indices, :w) || throw(ArgumentError("the ParticleList has no variable w"))
+    return ChargeWeights{T, typeof(pl)}(pl, only(pl.indices.w))
+end
+
 function ParticleList(h5::H5DataStore, path::AbstractString = "/")
     group = h5[path]
 
@@ -151,7 +239,15 @@ function ParticleList(h5::H5DataStore, path::AbstractString = "/")
 
     params = NamedTuple{paramkeys}(paramvals)
 
-    ParticleList(particles; variables = variables, parameters = params)
+    skwargs = if haskey(group, "species")
+        s = group["species"]
+        (species = Species.(Symbol.(read(s["name"])), read(s["charge"]), read(s["mass"])),
+            speciesindex = read(s["index"]))
+    else
+        (;)
+    end
+
+    ParticleList(particles; variables = variables, parameters = params, skwargs...)
 end
 
 function ParticleList(fpath::AbstractString, path::AbstractString = "/")
@@ -175,6 +271,12 @@ function h5save(h5::H5DataStore, p::ParticleList; path::AbstractString = "/")
         params[string(key)] = p.params[key]
     end
 
+    sgroup = _create_group(group, "species")
+    sgroup["name"] = [String(s.name) for s in p.species]
+    sgroup["charge"] = charge.(p.species)
+    sgroup["mass"] = mass.(p.species)
+    sgroup["index"] = p.speciesindex
+
     group["particles"] = p.list
 end
 
diff --git a/src/species.jl b/src/species.jl
new file mode 100644
index 0000000..fc454d2
--- /dev/null
+++ b/src/species.jl
@@ -0,0 +1,54 @@
+"""
+    Species(name::Symbol, charge, mass)
+    Species{T}(name::Symbol, charge, mass)
+
+A particle species: a name, a charge `q` and a mass `m`, shared by every particle of that species.
+`Species(name, charge, mass)` promotes the charge and the mass to their common type, and
+`Species{T}(s::Species)` converts an existing species to the element type `T` of a
+[`ParticleList`](@ref).
+
+The mass must be positive: a mass that is not, `NaN` included, throws an `ArgumentError`, and the
+conversion to `T` checks again, so a mass that rounds to zero in `T` throws too. The charge may be
+zero.
+
+# Example
+
+```julia
+s = Species(:e, -1.0, 0.5)
+s.charge   # -1.0
+s.mass     # 0.5
+```
+"""
+struct Species{T}
+    name::Symbol
+    charge::T
+    mass::T
+
+    function Species{T}(name::Symbol, charge, mass) where {T}
+        m = convert(T, mass)
+        m > zero(T) ||
+            throw(ArgumentError("the mass of species $name must be positive, got $m"))
+        new{T}(name, convert(T, charge), m)
+    end
+end
+
+function Species(name::Symbol, charge, mass)
+    Species{promote_type(typeof(charge), typeof(mass))}(
+        name, charge, mass)
+end
+
+Species{T}(s::Species) where {T} = Species{T}(s.name, s.charge, s.mass)
+
+"""
+    charge(s::Species)
+
+Return the charge `q` of the species `s`.
+"""
+charge(s::Species) = s.charge
+
+"""
+    mass(s::Species)
+
+Return the mass `m` of the species `s`.
+"""
+mass(s::Species) = s.mass
diff --git a/test/particle_list.jl b/test/particle_list.jl
index 4311cd7..96eed2f 100644
--- a/test/particle_list.jl
+++ b/test/particle_list.jl
@@ -101,3 +101,224 @@ plh5 = ParticleList(h5file)
 @test plh5[1, 1] == pl[1, 1]
 @test plh5[1, :] == pl[1, :]
 @test plh5[:, 1] == pl[:, 1]
+
+# A function barrier with concrete arguments, for the allocation assertions: the first call warms
+# the method and the second is measured.
+allocations(f::F, a::A) where {F, A} = (f(a); @allocated f(a))
+allocations(f::F, a::A, i::I) where {F, A, I} = (f(a, i); @allocated f(a, i))
+
+_allocating(pl, i) = zeros(i)
+_chargeweight(pl, i) = chargeweights(pl)[1, i]
+
+# the species table, the per-particle index, and the accessors
+@testset "species table and accessors ($T)" for T in (Float32, Float64)
+    np = 5
+    nd = 3
+    vars = (x = 1, v = 2, w = 3)
+    z = rand(T, nd, np)
+    specieslist = [Species(:e, -1, T(0.5)), Species(:i, 2, 4)]
+    speciesindex = [1, 2, 1, 2, 2]
+    pl = ParticleList(z; variables = vars, species = specieslist, speciesindex = speciesindex)
+
+    # the table is converted to the element type of the list, and the index stays `Int`
+    @test pl.species isa Vector{Species{T}}
+    @test pl.speciesindex isa Vector{Int}
+    @test pl.speciesindex == speciesindex
+
+    # every particle reads its own charge, mass and species
+    @test [charge(pl, i) for i in 1:np] == T[-1, 2, -1, 2, 2]
+    @test [mass(pl, i) for i in 1:np] == T[T(0.5), 4, T(0.5), 4, 4]
+    @test [species(pl, i).name for i in 1:np] == [:e, :i, :e, :i, :i]
+    @test charge(pl, 1) isa T
+    @test mass(pl, 1) isa T
+    @test species(pl, 1) isa Species{T}
+
+    # the accessors infer, and allocate nothing after a warm call
+    @test @inferred(charge(pl, 1)) == T(-1)
+    @test @inferred(mass(pl, 1)) == T(0.5)
+    @test @inferred(species(pl, 1)) == Species(:e, -1, T(0.5))
+    @test allocations(charge, pl, 1) == 0
+    @test allocations(mass, pl, 1) == 0
+    @test allocations(species, pl, 1) == 0
+    @test allocations(_allocating, pl, 1) > 0
+
+    # an `Int32` index works, on the accessors and in the keyword
+    @test charge(pl, Int32(2)) == charge(pl, 2)
+    @test mass(pl, Int32(2)) == mass(pl, 2)
+    @test species(pl, Int32(2)) === species(pl, 2)
+    @test ParticleList(z; variables = vars, species = specieslist,
+        speciesindex = Int32[1, 2, 1, 2, 2]) == pl
+
+    # the list copies the index, so a later write to the caller's vector leaves it valid
+    idx = copy(speciesindex)
+    pli = ParticleList(z; variables = vars, species = specieslist, speciesindex = idx)
+    idx[1] = 7
+    @test pli.speciesindex == speciesindex
+
+    # construction throws in three cases; the range is checked on both sides
+    @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist,
+        speciesindex = [1, 2, 1, 2])
+    @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist,
+        speciesindex = [1, 2, 3, 1, 2])
+    @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist,
+        speciesindex = [1, 2, 1, 2, 0])
+    @test_throws ArgumentError ParticleList(z; variables = vars, species = specieslist)
+
+    # equality compares the species table and the index
+    @test ParticleList(z; variables = vars, species = specieslist,
+        speciesindex = speciesindex) == pl
+    @test ParticleList(z; variables = vars,
+        species = [Species(:e, -1, T(0.5)), Species(:i, 3, 4)],
+        speciesindex = speciesindex) != pl
+    @test ParticleList(z; variables = vars, species = specieslist,
+        speciesindex = [2, 1, 1, 2, 2]) != pl
+
+    # a list built without a species holds one species, with charge and mass one
+    pld = ParticleList(z; variables = vars)
+    @test pld.species == [Species(:default, one(T), one(T))]
+    @test pld.speciesindex == ones(Int, np)
+    @test charge(pld, 1) == one(T)
+    @test mass(pld, 1) == one(T)
+    @test species(pld, 1).name === :default
+    @test vec(collect(chargeweights(pld))) == collect(pld.w)
+
+    # an empty species table with np ≥ 1 throws through the out-of-range check
+    @test_throws ArgumentError ParticleList(z; variables = vars, species = Species[])
+
+    # an empty list builds, with and without a species
+    z0 = rand(T, nd, 0)
+    @test length(ParticleList(z0; variables = vars, species = Species[])) == 0
+    @test length(ParticleList(z0; variables = vars)) == 0
+
+    # the conversion to the list's element type re-checks the mass
+    if T === Float32
+        @test_throws ArgumentError ParticleList(z; variables = vars,
+            species = [Species(:n, 1.0, 1e-50)],
+            speciesindex = ones(Int, np))
+    end
+end
+
+# the charge-weighted weights, and a deposit that uses them
+@testset "chargeweights" begin
+    np = 5
+    x = [0.4, 0.6, 1.4, 1.6, 1.9]
+    v = zeros(np)
+    w = [0.25, 0.5, 1.0, 2.0, 4.0]
+    specieslist = [Species(:e, -1, 0.5), Species(:i, 2, 4)]
+    speciesindex = [1, 2, 1, 2, 2]
+    pl = ParticleList(x, v, w; species = specieslist, speciesindex = speciesindex)
+
+    cw = chargeweights(pl)
+    @test cw isa AbstractMatrix{Float64}
+    @test size(cw) == (1, np)
+    # q_i · w_i by hand: (−1)(0.25), (2)(0.5), (−1)(1.0), (2)(2.0), (2)(4.0)
+    @test collect(cw) == [-0.25 1.0 -1.0 4.0 8.0]
+    @test @inferred(chargeweights(pl)[1, 1]) == -0.25
+    @test cw[1, Int32(2)] == 1.0
+    @test allocations(_chargeweight, pl, 1) == 0
+    @test_throws BoundsError cw[2, 1]
+
+    # a scalar `w` row and a one-row range agree, on one matrix
+    z = vcat(x', v', w')
+    pl_int = ParticleList(z; variables = (x = 1, v = 2, w = 3), species = specieslist,
+        speciesindex = speciesindex)
+    pl_rng = ParticleList(z; variables = (x = 1, v = 2, w = 3:3), species = specieslist,
+        speciesindex = speciesindex)
+    @test collect(chargeweights(pl_int)) ==
+          collect(chargeweights(pl_rng)) == [-0.25 1.0 -1.0 4.0 8.0]
+
+    # `w` spanning two rows, and a list without `w`
+    z4 = vcat(x', v', w', w')
+    @test_throws ArgumentError chargeweights(ParticleList(z4;
+        variables = (x = 1, v = 2, w = 3:4), species = specieslist, speciesindex = speciesindex))
+    @test_throws ArgumentError chargeweights(ParticleList(z; variables = (x = 1, v = 2),
+        species = specieslist, speciesindex = speciesindex))
+
+    # the deposit weights each particle by q·w, onto its nearest grid point
+    grid = [0.5, 1.5]
+    deposit = zeros(2)
+    for i in 1:np
+        j = argmin(abs.(grid .- pl.list[1, i]))
+        deposit[j] += cw[1, i]
+    end
+    # cell 1: (−1)(0.25) + (2)(0.5) = 0.75;  cell 2: (−1)(1.0) + (2)(2.0) + (2)(4.0) = 11.0
+    @test deposit == [0.75, 11.0]
+end
+
+# the species table and the index survive the HDF5 round trip
+@testset "HDF5 round trip of the species table ($T)" for T in (Float32, Float64)
+    np = 5
+    nd = 3
+    vars = (x = 1, v = 2, w = 3)
+    z = rand(T, nd, np)
+    specieslist = [Species(:e, -1, T(0.5)), Species(:i, 2, 4)]
+    speciesindex = [1, 2, 1, 2, 2]
+    pl = ParticleList(z; variables = vars, species = specieslist, speciesindex = speciesindex)
+
+    mktempdir() do dir
+        f = joinpath(dir, "species.h5")
+        h5open(f, "w") do file
+            ParticleMethods.h5save(file, pl)
+        end
+        @test ParticleList(f) == pl
+
+        # the datasets are written in the element type of the list, `index` as `Int`
+        h5open(f, "r") do file
+            @test eltype(read(file["species/charge"])) === T
+            @test eltype(read(file["species/mass"])) === T
+            @test eltype(read(file["species/index"])) === Int
+            @test read(file["species/name"]) isa Vector{String}
+        end
+    end
+end
+
+# a file whose species group is gone loads with the default species
+@testset "a file without a species group loads with the default" begin
+    np = 5
+    z = rand(3, np)
+    vars = (x = 1, v = 2, w = 3)
+    pl = ParticleList(
+        z; variables = vars, species = [Species(:e, -1, 0.5), Species(:i, 2, 4)],
+        speciesindex = [1, 2, 1, 2, 2])
+
+    mktempdir() do dir
+        f = joinpath(dir, "stripped.h5")
+        h5open(f, "w") do file
+            ParticleMethods.h5save(file, pl)
+        end
+        h5open(f, "r+") do file
+            delete_object(file, "species")
+        end
+        @test ParticleList(f) == ParticleList(z; variables = vars)
+    end
+end
+
+# the species group is written under `path`, beside `particles`
+@testset "the species group follows path" begin
+    np = 5
+    z = rand(3, np)
+    vars = (x = 1, v = 2, w = 3)
+    plA = ParticleList(
+        z; variables = vars, species = [Species(:e, -1, 0.5), Species(:i, 2, 4)],
+        speciesindex = [1, 2, 1, 2, 2])
+    plB = ParticleList(
+        z; variables = vars, species = [Species(:p, 3, 7), Species(:q, -4, 2)],
+        speciesindex = [2, 1, 2, 1, 2])
+
+    mktempdir() do dir
+        f = joinpath(dir, "nested.h5")
+        h5open(f, "w") do file
+            ParticleMethods.h5save(file, plA; path = "a")
+            ParticleMethods.h5save(file, plB; path = "b")
+        end
+        @test ParticleList(f, "a") == plA
+        @test ParticleList(f, "b") == plB
+
+        # a nested species group stripped from one path leaves the other alone
+        h5open(f, "r+") do file
+            delete_object(file, "a/species")
+        end
+        @test ParticleList(f, "a") == ParticleList(z; variables = vars)
+        @test ParticleList(f, "b") == plB
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index ce6c592..8e894f9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ const GROUPS = isempty(ARGS) ? ["core", "slow"] : ARGS
 if "core" in GROUPS
     @safetestset "Aqua" include("quality/aqua.jl")
     @safetestset "Particles" include("particle.jl")
+    @safetestset "Species" include("species.jl")
     @safetestset "Particle Lists" include("particle_list.jl")
     @safetestset "Defects" include("integration/defects.jl")
 end
diff --git a/test/species.jl b/test/species.jl
new file mode 100644
index 0000000..f6651ef
--- /dev/null
+++ b/test/species.jl
@@ -0,0 +1,39 @@
+using ParticleMethods
+using Test
+
+# a species carries a name, a charge and a mass, and the mass must be positive
+@testset "Species rejects a non-positive mass" begin
+    @test_throws ArgumentError Species(:n, 0.0, 0.0)
+    @test_throws ArgumentError Species(:n, 1.0, -1.0)
+    @test_throws ArgumentError Species(:n, 1.0, 0.0)
+    @test_throws ArgumentError Species(:n, 1.0, NaN)
+
+    s = Species(:n, 0.0, 1.0)
+    @test s isa Species{Float64}
+    @test s.name === :n
+    @test s.charge == 0.0
+    @test s.mass == 1.0
+end
+
+# the constructor promotes the charge and the mass, and the accessors read the fields
+@testset "Species promotes its fields" begin
+    e = Species(:e, -1, 0.5)
+    @test e isa Species{Float64}
+    @test charge(e) == -1.0
+    @test mass(e) == 0.5
+
+    i = Species(:i, 2, 4)
+    @test i isa Species{Int}
+    @test charge(i) === 2
+    @test mass(i) === 4
+end
+
+# Species{T}(s) converts an existing species, and re-checks the mass in T
+@testset "Species converts to another element type" begin
+    e = Species(:e, -1, 0.5)
+    @test Species{Float32}(e) isa Species{Float32}
+    @test charge(Species{Float32}(e)) === -1.0f0
+    @test mass(Species{Float32}(e)) === 0.5f0
+
+    @test_throws ArgumentError Species{Float32}(Species(:n, 1.0, 1e-50))
+end